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ABSTRACT 


LEYBOLD, HERBERT ARTHUR. Simulation of Service . ad Dis- 
tributions Utilizing Stationary, Non-Gaussian Random Time 

Histories. (Under the direction of FRANKLIN DELANO HART). 

# 

The objective of the study is to determine the feasi- 
bility of generating a non-Gaussian, stationary, random 
load time history which can be used to simulate service 
load histories during fatigue testing of structural parts. 

To achieve this objective, five random load time his- 
tories of arbitrary amplitude and each having a different, 
mathematically describable, non-Gaussian, amplitude prob- 
ability density function, were generated with the aid of a 
CDC 6000 series digital computer and FORTRAN IV programming. 
Cumulative peak distributions were obtained from each of 
the resulting random load time histories and compared with 
cumulative peak service load distributions. Five random 
load time histories are Investigated in order to assure 
that the best possible duplication of the cumulative peak 
service load distribution is obtained. 

Of primary importance In the Investigation was the 
matching of predicted, that is, digitally generated and 
service load cumulative peak distributions. Of secondary 
importance was the technique used for digitally generating' 
the five different, non-Gaussian, random load time 
histories. 

Cumulative peak distributions of two of the generated 
random time histories approximated the cumulative peak 



service distributions for aircraft . * The remaining three 
cumulative peak distributions approximated the cumulative 
peak service distributions of ground transportation, that 
is, car, truck, rail, overhead travelling crane, etc. 
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1 . INTRODUCTION 


Over the past several years the trend has been’ jtowards 
random load fatigue testing of aircraft structural parts 
rather than the constant amplitude and programmed ccaistant 
amplitude fatigue testing which has been going on for /years . 
The advent of more sophisticated fatigue testing machines 
capable of applying complex load-time histories which more 
nearly simulate actual service experience has been th% major 
factor towards this trend. Within this framework of random 
load fatigue testing there has evolved two schools of 
thought on simulating aircraft service experience. The 
first utilizes a quasi-stationary Gaussian random process to 
simulate service experience. This is accomplished by pro- 
gramming, in a random fashion, the root-mean-square level 
of a Gaussian random process. The second utilizes a 
stationary non-Gaussian random process to simulate service 
experience. It is the purpose of this thesis to show that 
this second approach is feasible. > — - 

Toward this end a random number generator was used to 
generate random numbers having a uniform or constant ampli- 
tude probability density function between the limits of 
zero and one. These numbers were then shaped such that the 
amplitude probability density of a discrete random time 
history composed of these random numbers had one of the 
following density functions: Poisson, Binomial, Log-Formal, 

Weibull, and Exponential. These density functions were 
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arbitrarily selected because they are widely discussed in 

the literature and are mathematically describable. The 

* 4 ' • . . 

'i resulting random time histories sire analyzed to determine 
their peak statistics. It is generally believed that the 
peak statistic is the significant fatigue inducing parameter 
of a random time history. It is also the statistic most 
frequently obtained during service experience. Cumulative 

j* 

peak distributions are computed and compared with cumulative 
peak service distributions. For aircraft the distribution 
is a straight line on semi-log paper. 

It is shown that the random time histories having 
Exponential and Log-Normal amplitude probability density 
functions can be used to simulate the service load experi- 
ence of aircraft. It is also shown that the other three 

~ * ■ 

random time histories having Binomial, Poisson and Weibull 
amplitude probability density functions can be used to 
simulate service load experience of ground transportation 
„such as cars, trucks, railroads, and overhead travelling 
cranes * 
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2. REVIEW OP LITERATURE 

It. ha3 become increasingly more apparent' in recent 

( ■ j? 

years that more and more fatigue tests are being conducted 

> i 

using random loads. There appears to be several reasons 
for thisftrend. First, considerably more research has been 
done in recent years in this area as evidenced by the number 
of articles being published. These articles have undoubtedly 
influenced the fatigue engineers philosophy or outlook on 

t Jjh 

fatigue Resting. Secondly, the fatigue engineer has come 
to realize that service loadings are complex loadings and 
that he ’is incapable of predicting fatigue life under their 
influence. Hence he must resort to a much more accurate 
simulation of the service loadings in order to determine 
fatigue life. Finally, he now has at his disposal testing 
equipment- which will permit him to apply any complex load 
history he desires to test specimens and structures. The 
above mentioned items have brought us to the current state 

’* -S' 

of the art in random load fatigue testing. Random process 
theory upon which random load fatigue testing is based cam 
be found^ln such texts as Bendat and Piersol (I 966 ), 

Crandall and Mark (1963), Davenport and Root (1958) and 

v';’- 

Robson (I 96 I 4 .). 

•s 

Swaifson (1968) presents an excellent review of the 
literature on random load fatigue testing. In this refer- 
ence h-e. starts with the testing- to develop the well known 
S-N curve and goes on to the block tests wherein constant 
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amplitude cycles are stratified at various stress levels in 
an arbitrary sequence so as to approximate the relative 
frequency of service*, loads from the'service load spectrum. 

In addition he discusses the actual random process and/or 
an analogous random .process to the one encountered in 
service which may be either stationary, quasi- stationary or 
non- stationary and result in the same load frequency as 
given by the service load spectrum. Swanson also Indicates 
that there are basically two approaches to fatigue testing 
under random loading: the "digital” approach and the 

analogue approach. /The digital approach will necessarily 
result in a random process that is stationary in nature 
while the analogue approach may result in either a sta- 
tionary or quasi- st at ionary process. 

The author (1963) was involved in obtaining statistics 

for the conduct of bpth block and random type tests. Based 

* 

on these statistics-^ Leybold and Naumann (1963) conducted 

tests to evaluate the effects of both block and random type 

•***•*.- * ■ 

tests on the fatigue life of aluminum alloys. Leybold ( 1963 ) 

- k 
v. 

noted that the primly statistic used to conduct fatigue 
tests was the peak statistic. A mathematical method for 
predicting the. peak distribution or peak statistic has been 
developed by Rice (1954) for a Gaussian, stationary random 

V 

process. Broch (1963 & 1968) discussed the Work developed 
by Rice and applies ^.t to the case of f atigue under random 
loading. Schijve' (1963) tried to determine which statistic 
of a random process was most significant with regard to 
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fatigue. In most cases the random process investigated and 
used for testing was Gaussia^ and stationary in nature. Two 
primary reasons for using this type of a random process 
were: (1) a Gaussian stationary random process is easily 

generated with a commerciall^available random noise gener- 

"* fc 

ator and (2) the peak distribution of a Gaussian, stationary 

4 . f 

random process was easily commuted based on Rice's work. 

£ 

However, the drawback to uslaft this type of a random process 

is that service load data sho^r that they are rarely station- 

&' 

ary in nature and may or may, 'not be Gaussian in nature. 

' * & 

Based on. this fact, Swanson, (1963 & 1965) has proposed the 

• 4 . 

Gaussian non- at at ionary r.andw process for fatigue testing. 

Using this method of testing ^the rms of a stationary random 

- ej;. 

process is programmed to vary' in a random fashion such that 
the resulting peak distribution will match that of service 
load data. Service load dat^are usually presented in the 
form of cumulative peak distributions, that is, cumulative 
number of peaks exceeding a gfyen level. Service load 
cumulative peak distributions , for aircraft usually plot as 

V. \ 

straight lines on semi-log pgper. To cite just two of the 

-f' , 

publications available with aircraft service load data we 

- -V \ 

have Coleman (1968) who presents his data in terms of g's 
(acceleration due to gravity and Bullen (1967) who pre- 
sents his data in terms of gust velocities. 

The quasi-st at ionary - random fatigue test proposed by 
Swanson is one way of simulating the. service load peak 
distribution. Another method also proposed by Swanson (1968) 
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was to find a stationary non-Gaussian, random process, which 
exhibits a cumulative peak distribution the '» same as an air- 

A 

> V > - f 

craft service load cumulative peak distribution, that is, a 

■ V . 

straight line on semi-log paper. Thus two! different ways 
have been proposed for simulating service ^dad time 
histories: (1) the .quasi -stationary methdd^where the 

random process is Gaussian in nature and the root-mean- 
square levels are programmed in a random fashion and 
(2) the stationary method, that is, constant root -mean- 
square level, where the random process is non-Gaussian in 
nature. This thesis deals with the latter, method. 
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3. GENERATION AND ANALYSIS OP RANDOM TIME HISTORIES 

•J 

3,1 General . „ 

In this chapter random numbers are generated and,, 
shaped so that the resulting discrete random time hi Story 

has a non-Gaussian amplitude probability density function. 

\ • . 

Actual and theoretical amplitude probability density func- 
tions are compared. In addition peak and cumulative peak 

i p. 

distributions are calculated. The computations are broken 
down into three parts as shown in the next three sections. 
The computer programs for carrying out these computations 

-i »*» 

can be found in the Appendices, sections 8.1 and 8.2. More 
will be said about these programs in Chapter I 4 . on computer 
programs . 

Five different amplitude probability density- functions 
were used in order to insure that the best possible agree- 
ment between the digitally generated data and the service 
load data would be attained. 

It should be noted that the data presented are only a 
small portion of the data generated with the ai 1 of. 
computer. Parametric studies were made using all five*/ 
amplitude probability density functions. The resulting^' 

data were just too voluminous to include in this thesis. 

> , 

Only those data which show trends or agreement with service 
load data have been included. 
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3.2 Generation of Random Numbers 
Random numbers bounded by the limits zero and one were 
generated with' a standard random number generator having a 
uniform "br constant amplitude probability density function. 
The numbers were then shaped such that the resulting dis- 

i ^ 

Crete ra'iidom time history would have some arbitrarily 
selected -amplitude probability density function. The 
followin'^ amplitude probability density functions were 
chosen because they are mathematically describable and 
discussed thoroughly in the literature: Poisson, Binomial, 

Log-Norm|:l, Weibull and Exponential. The procedure for 
digital]^- generating the shaped random numbers was as 
followsjf- Let F(x) be the amplitude probability density 
function and CF(x) the cumulative amplitude probability 
density function. CF(x) is related to the amplitude 
probability density function by the relation 

dCF(i) = F(x)dx (3.1) 

or by integrating 

CF(x) = | F(a)da (3.2) 

U — 00 

which is bounded by the limits zero and one when x takes 
on the'. values of and +•, respectively. Thus, random 

numbers having a uniform probability density distribution 
with limits between zero and one can be generated with a 
random number generator and substituted for CF(x) in 
equation (3.2) . This equation can then be evaluated . 
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for x. As many random numbers can be generated in this 
.manner as are needed and they will necessarily have the 
.desired amplitude, probability density function. For this 

work each random time history consisted of 10000 random 

• * 

numbers. 

The five random time histories generated in this 
manner had the following amplitude probability density 
functions: Poisson, Binomial, Log-Normal, Weibull and 

Exponential. In the first of these time histories, whose 
amplitude probability density function is discrete, the 
following recurrence formula was used to compute the 
cumulative probability density: 

F(x + 1) =7^TPU) (3.3) 

The cumulative frequency, CF(x), was compared with a 
digitally generated uniformly distributed random number Y 
having a value between zero and one. 

The random variable x was found by determining 
when Y fell between the limits CF(x) and CF(x + 1), 
that is, 

• " c x 

. v Y = X p (*) (3-k) 

The second timjfiL history wa? generated in a similar 
manner using the binomial recurrence formula: 

* ’VerV 

F(x + 1) ■ ~ - p (*J 

(x + 1) q 


(3.5) 
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The third time history or log-normal one was generated 
by solving equation (3.2) with a numerical approximation of 
the integral (see section 8. 2. 2,3). 

In the case of the last two time histories equa- 
tion (3.2) could be integrated directly and resulted in 
the following random variables ’ xft For the Weibull prob- 
ability density function: 

F(x) = apx a " 1 exp(-Px a ) (3.6) 

and f 

CF(x) = 1 - exp(-Px a ) (3.7) 


solving equation (3.7) for x we get 

x = log e [l ^CF(x)]j a 

For the Exponential probability density function: 

F(x) =0expt-ex) (3.8) 

* 

and 

CF(x) = 1 - eip(-0x) (3.9) 

solving equation (3.9) for x we get 

x =-| log e [l - CF(x)] 

It should be noted that the Exponential distribution is a 
special case of the Weibull distribution when a.= 1. 
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The first 100 random numbers generated for each of the 
five different amplitude probability density functions used 
are shown in Figures 3.. 1-3.10. The mean, r pot-mean- square 


and standard deviation are shown for each time history. 

- i 

The mean indicates the average or expected value of the 
random variable, say x. The mean is mathematically 


expressed as 


E(x) 


p ® 

I xF(x)dx 



( 3 . 10 ) 


The root -me an- square is a measure of the amplitude of a 
time history and is the square root of the mean of the 
squares of a random variable x. The mean square is 
mathematically expressed as 


3{x 2 ) 




x 2 F(x)dx 


(3.11) 


The standard deviation is a measure of the dispersion from 
the mean and is the square root of the variance of a random 
variable x. The variance is defined mathematically as 

a 2 = E - E(x)] 2 } - v (3.12) 

Crandall and Mark (1963) show that this can be reduced to 

o 2 = E(x 2 ) - [E(x)] 2 - (3.13) 

" 

4 

All of the above mentioned properties, were calculated 

if '* 

\ for each of the five random time histories generated. The 
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Figure 3.1. Time history of random numbers having a Weibull 
amplitude probability density: a = 2j $ = 2 
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Figure 3.2, Time history of random numbers having a Weibull 
amplitude probability density; a =5* p = 5 
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Figure* 3. 3* Time history of random numbers having an 
Exponential amplitude probability density: 
8 = 1 
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Figure 3.I4.. Tim© history of random numbers having an 
Exponential amplitude probability density: 
0 = 2 
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Figure 3.5. Time history of random numbers having a Poisson 
amplitude probability density: p - 10 
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Figure 3.7. Time history of random numbers having a 
Binomial amplitude probability density: 
n = 10; p = 0.1 



13 



Time 

v 4 

.•Figure 3.8. Time history of random numbers having a 
Binomial amplitude probability density: 
n = 500; p = 0.1 
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Figure 3.10. 


Time history of randan numbers having a 
Log-Normal amplittde probability density? 

ef * 0.05 

log 
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results of these calculations are shown in Table 3.1. In 
order to check these values with the theoretical ones, the 
author scanned the literature to obtain the theoretical 
equations for the mean, mean-square and variance of the 
five density functions used. To his dismay all were not 
readily available. Thus, the basic statistical equations 
and density functions were used to derive the desired 
properties. These derivations are shown in Appendix 8.5 
and the results are summarized in Table 3 . 2 . The theo- 
retical properties were calculated and showed good agree- 
ment with the properties tabulated in Table 3 . 1 . 

3.3 Calculation and Comparison of Computed 
and Theoretical Histograms 
After each of the five random time histories Iwere 
digitally generated, the amplitude scale of each was 
divided up into equal intervals and the random numbers 
falling in each interval were counted for all lb, 000 num- 
bers of each time history. The resulting frequency of 
occurrence of the random numbers were plotted as histograms 
in Figures 3.11-3.22. In addition the theoretical histo- 
grams were also plotted and compared with the .actual 
histograms in these figures. It can be seen that good 
agreement was obtained. A chi-squared test was performed 
to compare the actual histogram with the theoretical one. 
In all cases the calculated chi-square value ^as well 

r y '■ 

below the theoretical chi-squared value at the. '95 percent 



Table 3.1. Comparison of computed and theoretical properties of density functions 





Table 3.2. : Summary of properties .of density functions 
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Figure 3.1ij.. Actual and theoretical histogram for an 
Exponential amplitude probability density: 






29 



Figure 3.15. Actual and theoretical histograms for a ,# 
Poisson amplitude probability density: 

p = 10 
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Figure 3.16. Actual and theoretical histograms for a 
Poisson amplitude probability density: 

\i = 66 
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l theoretical histograms for a 
amplitude probability density: 


0.1 
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Figure 3.19. 


Actual and theoretical hiatograma for a 
Log-Normal amplitude probability denaity 

®log * °* 1 
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Figure 3.20. Actual and theoretical histograms for a 
Log>Normal amplitude probability density: 

°fog s °*°* 
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Figure 3.21. Actual and theoretical logarithmic histograms 
for a Log-Nomal amplitude probability 

density: a? = 0.1 




Frequency of occurrence 
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Figure 3 . 22,6 Actual and theoretical logarithmic histograms 
for a Log-Normal amplitude probability 
density: =. 0,0$ 
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confidence level. This exercise was done in order to show 
that the digitally generated random time histories do 
indeed have the desired amplitude probability density 
functions. 

As indicated in the title of this thesis all time 
histories generated we-re to be non-Gaussian in nature. 
However, certain values for the various parameters in the 
density functions used tend to make the distributions 
approach a Gaussian distribution. It is well known that 
some of these distributions can be approximated with a 
normal or Gaussian distribution. This fact becomes obvious 
when looking at Figures 3.11-3.22. The a and (3 param- 
eters of the Weibull distribution were each varied in unit 
increments from 2 to 5. Figure 3.11 where a = 2, 0 = 2 
shows a highly non-Gaussian shape while Figure 3.12 with 
a = £, p = £ shows a shape approaching Gaussian. The 9 
parameter of the Exponential distribution was varied from 
l/!f through 2 by doubling each 0 in succession. This 
distribution will always be highly non-Gaussian as shown 
in Figures 3.13 and 3.14. The Poisson distribution tends 
toward the Gaussian distribution even for small mean 
values i* as illustrated in Figures 3.15 and 3.16. The 
Binomial distribution tends to be highly non-Gaussian for 
small n (see Figure 3.17) and Gaussian for large n 

(see Figure 3.18). In the Log-Normal distribution, the 

2 

parameters u, and a were varied as follows: 

• log v * log 
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2 

= 0, 1; CT i 0 g = 2, I* 0.5, 0.1, 0.05. A horizontal 

translation of the distribution occurred when ^i 0 g was 
increased from zero to one. Figures 3.19 and 3.20 show the 
highly non-Gaussian shape of the distribution while 
Figures 3.21 and 3.22 show the , Gaussian shape when 
plotted on semi-log coordinates. 

3 .I 4 . Determination of Peak and Cumulative Peak 
Distributions for Both Maximum and 
Minimum Peak Values 

As with the histograms, the amplitude scale of the 
time histories were divided up into equal Intervals and the 
peaks, both maximum and minimum, were counted in each 
interval. Maximum peaks and minimum peaks are defined as 
shown in Figure 3 . 23 . The maximum peak distributions are 
plotted in Figures 3* 21*.— 3 .33 while the minimum peak dis- 
tributions are shown in Figures 3 * 31*-- 3 • 14-3 ^ Since service 
load data are based on maximum peaks, we are primarily 
interested in .the maximum peak distributions of the gen- 
erated time histories. However, it is also of Interest to 
see whether or not the peaks, both maximum and minimum, 
are distributed symmetrically about the mean of a random 
time history* Thus the minimum' peak distributions are 
also plotted. One can check for symmetry by visualizing 
a maximum peak distribution overlaid on a minimum peak 
distribution. The two distributions will intersect at the 




if .0 



Class Mark 


Figure 3.2I4.. Maximum peak distribution of a time history 
having a Weibull amplitude probability 
density: a = 2; p = 2 
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.Fig-are 3.25. Maximum peak distribution of a time history 
% having a Weibull amplitude probability 

density; a = 5; P = 5 
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Figure 3.26, Maximum peak distribution of a time history 
having an Exponential amplitude probability 
density: 9 « 1 
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Figure 3.27. Maximum peak distribution of a time history 
having an Exponential amplitude probability 
density: 9=2 
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Figure 3.28. Maximum peak distribution of a time history 
having a Poisson amplitude probability 
density: = 10 
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Figure 3.3 0. Maximum peak distribution of a time history 
having a Binomial amplitude probability 
density: n = 10; p * 0.1 
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Figure 3.31. Maximum peak distribution of a time history 

having a Binomial amplitude probability ■ * * 
density: n - 500 ; p » 0.1 
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Figure 3.3 2. Maximum peak distribution of a time history 
having a Log-Normal amplitude probability 

density: a? ~ 0,1 

log 
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Figure 3.33* Maximum peak distribution of a time history 
having a Log-Normal amplitude probability 

density* a f C g ~ 0 , 0 $ 
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Figure 3*36. Minimum peak distribution of a time history 
having an Exponential amplitude probability 
density!; 0 = 1 
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Figure 3.37. Minimum peak distribution of a time history 
having an Exponential amplitude probability 
density: 0=2 
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Figure 3.36* Minimum peak distribution of a time history 
having a Poisson amplitude probability 
density? ji « 10 
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Figure 3.39. Minimum peak distribution of a time history 
having a Poisson amplitude probability 
density? u. ~ 66 






Frequency of oceurr 



Figure 3.40, Minimum peak distribution of a time history 
haying a Binomial amplitude probability 
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distribution of a time histo: 
omial amplitude probability 

= 500 ; p = 0.1 
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* * * 

Figure 3.1+2. Minimum peak distribution of a time history 
having a Log-Normal amplitude probability 

density: a fog = 
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Figure 3.l;3o Hinimian peak distribution of a time history 
having a Log-Normal amplitude probability 

density; ff i 0 g ® 
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histories presented in Figures 3.1-3.10. It can be seen 
from the peak distributions which are not symmetrical about 
the mean that there are a greater number of maximum peaks at 
the lower amplitudes "and a correspondingly smaller number 
of peaks at the higher amplitudes. This is exactly the 
same condition that Is encountered in aircraft. There are 
a relatively high number of gusts of low velocity with a 
corresponding decrease in the number of gusts at high 
velocity. Maximum peaks are measured as the number of 
exceedances of a given amplitude level above the mean or Ig 
stress of an aircraft. Thus the peaks measured in service 
are presented as a cumulative distribution. Thus in this 

investigation c\naulative peak distributions ..were computed 

.... " • ' 

for both maximum and minimum peaks and are presented in 

Figures and Figures 3«5 J 4~3®63» respectively. 

The maximum cumulative peak distributions are seen to start 
rolling off at the mean value of the time history. Only 
in the case of the Exponential and Log-Normal density 
function do these distributions fall off from the mean on 
a straight line. This fact is most significant in that 
this is the manner in which service data behaves. More 
will be said about this point in chapter 5. 
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Figure 3.U4* Log of cumulative maximum peaks of a time 
history having a Weibull amplitude 
probability density: a = 2; p = 2 
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Figure 3.45. Log of cumulative maximum peaks of a time 
history having a Weibull amplitude 
probability density? a = P = 5 
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Figure 3.4-6. Log of emulative maximum peaks of 
history having an Exponential ampli 
probability- density: 0 = 1 
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Figure 3.i;7. Log of cumulative maximum peeks of a time 
history having an Exponential amplitude 
probability density: 0-2 
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Figure 3 . 14 - 9 . Log of cumulative maximum peaks of a time 
history having a Poisson amplitude 
probability density: p = 66 
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Figure 3.j?l . Log of cumulative maximum peaks of a time 
history having a Binomial amplitude 
probability density: n - £00* P 35 0*1 
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Figure 3.52. Log of cumulative maximum peaks of a time 
history having a Log-Normal amplitude 

probability density: or| og “ 0.1 
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Figure 3.53* Log of cumulative maximum peaks of a time 
history having a Log-Normal amplitude 

probability density: cr? = 0.05 
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Figure 3. £5* I<og of cumulative minimum peaks of a time 
history having a Weibull amplitude 
probability density: a = p = 5 
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Figure 3.57* Log of cumulative minimum peaks of a time 
history having an Exponential amplitude 
probability density; 9-2 
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Figure 3.58. Log of cumulative minimum peaks of a time 
history having a Poisson amplitude 
probability density: p = 10 
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Figure 3.59. Log of cumulative minimum peaks of a time 
history having a Poisson amplitude 
probability density: y. = 66 
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Figure 3* 60< Log of emulative minimum peaks of a time 
history having a Binomial amplitude 
probability density: n * 10 j p = 0.1 
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Figure 3*61. Log of cumulative minimum peaks of a time 
history having a Binomial amplitude 
probability density? n = 500 ; p * 0.1 
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Figure 3*62. Log of emulative minimum peaks of a time 
history having a Log-Normal amplitude 
probability density! a 2 =0.1 
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Figure 3.63. Log of emulative minimum peaks of a time 
history having a Log-Normal amplitude 

probability density: of =0.05 

log 
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14. . COMPUTER PROGRAMS 
v I4..I General 

The "programs presented here were developed In order to 
digitally generate stationary, non-Gaussian random time 
histories which have mathematically describable amplitude 
probability density functions. The programs were written 
in FORTRAN IV language and ran on a CDC 6000 series digital 
computer. Compile time for the programs was something less 
than 10 seconds with an execution or run time of approxi- 
mately 30 seconds per each parameter variation, 

I4..2 Program Listings and Typical Block Diagram 

The main programs are listed in Appendix 8.1 as 
follows: 

8.1.1 Binomial Distribution 

8.1.2 Poisson Distribution 

8.1.3 Weibull Distribution 

8.1.U- Exponential Distribution 

8.1,5 Log-Normal Distribution 

Subroutines are listed in Appendix 8.2 as follows: 

8. 2. 1.1 Random Number Generator RANF 

8. 2. 1.2 Plotting Subroutines DDIPLT 

8. 2/2.1 Function XWEIB 

8.2.2. 2 Function XEXP 

8 . 2. 2. 3 Function XLOGNL 

A typical block diagram for these computer programs is 
shown in Appendix 8.3. 
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U.3 Input Data 

Required input data are defined in comment cards near 
the beginning of each program listing (see section 8.1) . 

I 4 ..I 1 - Output Data 

The output of the programs which have been discussed 
in section 3 are summarized in plot form. Specifically, 
the outputs are time histories, histograms, peak distribu- 
tions and cumulative peak distributions. 



5. COMPARISON OP CUMULATIVE PEAK DISTRIBUTIONS 
WITH SERVICE LOAD DATA 

5.1 Applicability 

Thus far the discussion has been centered around the 
simulation of service load data for aircraft. In most 
cases this type of data is presented for gust loads on 
either a ”g" (acceleration due to gravity) or velocity 
basis. The data are usually plotted on semi-log paper with 
the cumulative number of exceedances being plotted on the 
log scale and the "g” or velocity level on the linear 
scale. In first approximation this curve plots as a 
straight line. This is illustrated in Figure 5.1 which has 
been taken from Bullen (1967). 'It should be pointed out 
that service load data are based on maximum peak values. 

If this typical gust load spectrum is compared with the 
digitally generated cumulative peak distributions, it can 
be seen that the data presented In .Figures 3.4&* '3.47* 3-52 
and 3.53 could very easily be used to approximate the 
service load data by appropriate scaling of the abscissa. 

As an example, the digitally generated data of Figure 3*52 
has been replotted in Figure 5*2 along with the gust load 
data of figure 5.1. The abscissa of Figure 3.52 was scaled 
to agree with the abscissa of Figure 5.1. It can be seen 
that the digitally generated data agrees very well with the 
service load data. In the case of the gust load data pre- 
sented by Coleman (1968), where the cumulative frequency of 
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Figure £.2. Comparison of digitally generated data and 
service load data 
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exceeding an acceleration of a given level is given on a 
per mile or per flight basis, all that is needed is an 
estimate of the total number of miles or flights to be 
flown in the aircraft lifetime. With this information 
these curves become comparable with other cumulative peak 
distributions found in the literature. 

With regard to the other cumulative peak distributions 
which were digitally generated, they may be compared favor- 
ably with service load data for other than aircraft. 

Prothro (I 968 ) and Jacoby (I 967 ) present service load data 
for ground transportation items such as cars, trucks, and 
railroads and also for traveling overhead cranes. 

Figure 5.3 shows typical loading spectra for motor 
cars where the load time history was measured over varying 
distances from 1 km to 1250 km. The shorter the measuring 
distance the more apt one is to encounter only one specific 
road condition, that is, a concrete highway, a tar surface, 
a country road, etc. One surface condition usually means 
the load is distributed normally and hence the cumulative 
distribution is rounded as shown by the two lower curves 
in Figure 5.3. The longer the measuring distance is the 
more one is apt to get a mixture of surface conditions in 
the measured interval and hence the distribution changes 
from a normal one to a straight line distribution as indi- 
cated by the upper two curves in Figure 5.3. It should be 
noted that the spring displacement at the intersection with 
the ordinate is the mean spring displacement, that is, the 
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Figure 



5.3* Typical motor car service loading spectrum 
(Jacoby 1967) 



displacement due to the weight of the car cn the spring®. 

It can be seen that these curves compare favorably with the 
cumulative maximum peak distributions shown in Figures 3 * Uk = 
3.53. 

Figure %.h is. a typical service load spectrum for the 
front axle of a truck, this time in terms of bending fhePieftt . 
Again, as was discussed in Figure 5.3, the distribution 
starts at the mean or static load position. The dished 
curve represents what happens if the loads on the front 
axle of the truck are such that the springs '’bottom out . ** 
This phenomenon is not observed with the digitally gen- 
erated distributions but the basic cumulative frequency 
curve is similar to some of those shown in Figures 
3.53. 

Figure 5.5 shows a typical load spectrum for a rail- 
road car travelling at 65 mph. The measured distance was 
173 miles. Although one might assume that the surfaoe con- 
dition of a track would be constant, Prothro (I960) points 
out that track is classified as smooth and rough. In fact 
7 I 4 . miles of the 173 measured miles was classified as 
smooth. Since we can classify track as to rough, smooth 
or somewhere in between it is logical to assume that in 
the 173 measured miles there was a mixing of the various 
track conditions. Using the same reasoning as for motor 
cars we would expect a straight line distribution. We do 
indeed get a straight line distribution above the mean 
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Figure 5-U-- Typical truck front axle service loading 
spectrum (Jacoby 1967) 
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stress which is approximately 0.78g. Again this curve com- 
pares favorably with some of the digitally generated cumu- 
lative distributions shown in Figures 3.UU-3.53. 

Finally, Figure 5.6 shows a typical service load 
spectrum for a travelling overhead crane girder* Those 
different loadings are superimposed on this curve, all 
starting out from the mean stress. The basic loading forms 
the familiar curved distribution Indicating a somewhat 
Gaussian distribution of loads. This should have been 
anticipated since the travelling surface of the crane would 
not be expected to vary much in roughness and the measured 
distance is small. The vibrations Indicated are of small 
magnitude and relatively high frequency while the overload 
stresses are of larger magnitude and account for relatively 
few loadings. As with the other specimens discussed the 
shape of the basic loading curve is similar to some of 
these shown in Figures 3.W4--3.53. 

At this point it is clearly obvious that service load 
data from many sources, aircraft, car, truck, railroad and 
crane can be simulated using digitally generated non- 
Gausaian stationary random time histories. Once the basic 
shape of the cumulative distribution is obtained the ordi- 
nates and abscissas of these distributions can be adjusted 
to coincide with the appropriate service load data. 

The comparison of the cumulative peak distributions of 
time histories having Wei bull, Poisson and Binomial 






probability density functions was not made with the service 
load data from ground transportation because of the varied 
number of cumulative peaks exhibited (see Figures 5»3-5«6). 
The numb er of peaks is a direct function of the sample 
length of the time history. In the present investigation 
the sample length was fixed arbitrarily at 10,000 random 
numbers which resulted in 3*336 maximum peaks. The sample 
length .could just as easily be adjusted to result in 5,000 
or 10,000 peaks or any other number of peaks so as to match 
the service load data. The primary purpose here was to 
show that the shape of the cumulative peak distribution 
matched the service load data. 

5.2 Utilization 

It has been pointed out previously that two methods 
have been proposed for simulating. In a random fashion, the 
service load histories of aircraft. The first method pro- 
poses the generation of a "quasi-st at ionary” random time 
history whose emulative peak distribution simulates that 
found in service. The second method, described herein, 
proposes a stationary, non-Gaussian random time history 
whose cumulative peak distribution also simulates that 
found in service. Either or both methods can be used to 
conduct random load fatigue tests under simulated service 
conditions. 

Proponents of the quasi -stationary method of testing 
argue that the root -me an- square stress or load encountered 
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in service is varying continuously or is at least only 
constant for short periods of time. The use of a station- 
ary Gaussian time history from a random noise generator 
whose root-mean-square level is varied in a random fashion 
is employed to simulate the service conditions. The number 
of root-mean-square levels chosen to represent the condi- 
tions encountered ir. service is of primary importance. Too 
few levels will not yield the required service distribution. 
Too many will complicate the test setup. Of course, this 
method of testing necessarily uses an analogue signal as 
its input either to an electrodynamic or an electro- 
hydraulic shaker. It is also said that this method enables 
the wave form of the analogue signal to be preserved, that 
is, the frequency characteristics of the signal will be 
applied to the specimen as generated. However, Cl evens on 
and Steiner ( 1967 ) have shown that the shape of the power 
spectral during curve of a random time history, that is, 
the frequency content of the random signal, has little or 
no effect on the fatigue life. 

The stationary approach described herein, which is 
necessarily digital, can be cited for its ease of testing. 

A punched tape or punched card can be used to actuate a 
hydraulic servo valve to apply the peak loads in their 
proper sequence. The root-mean- square level is held 
constant so that no adjustment is required throughout a 
test. Wave form is lost with this method but as mentioned 
previously is considered of no consequence. 



9 c 

Analytically both methods result in the same emulative 
peak distribution. In addition, both methods apply the 
loads in a random failure. Based on a linear accumulation 
of damage both methods would necessarily result in the same 
fatigue life. 

It is the author's opinion that the digital or sta- 
tionary approach to testing is the more practical of the 
.two .proposed methods for the following reasons: (1) all 

the peak statistics are determined analytically prior to 
testing, (2) input to testing machine is in the form of 
punched or magnetic tape, and (3) a constant root-mean- 
square level is maintained throughout the test. However, 
it would be beneficial to conduct both stationary and 
quasi- stationary tests in order to provide additional 
incite into the important fatigue inducing parameters 
involved in fatigue testing under random loading. 
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6. CONCLUSIONS AND RECOMMENDATIONS 

Five random time histories having mathematically 
describable, non-Gaussian, amplitude probability density 
functions have been generated with the aid of a digital 
computer and their peak statistics determined. The fol- 
lowing conclusions have been drawn from an analysis of the 
data obtained: (1) Stationary random time histories having 

Exponential and/or Log-Normal amplitude probability density 
functions can be used to simulate the peak load distribu- 
tions encountered by aircraft in service. (2) Peak load 
distributions encountered in service by ground transpor- 
tation and handling equipment , that is, car, truck, rail 
and overhead crane, can be approximated using the remaining 
tliree time histories investigated, namely those having 
Weibull, Poisson, and Binomial amplitude probability 
density functions, (3) The resulting time histories can 
be used, to conduct stationary random and/or programmed 
block fatigue tests in which they simulate service load 
conditions. (I 4 .) Maximum and minimum peak distributions 
are unsymmetrical about the mean for non-Gaussian proba- 
bility density f met ions but approach a symmetrical con- 
dition as the probability density function approaches the 
normal or Gaussian distribution. 

It is recommended that a comparison based on fatigue 
testing be made between the proposed "stationary" method 
of simulating service conditions and the proposed 



"quasi-stationary" method of simulating the same conditions. 
These tests could provide additional incite into the 
important fatigue inducing parameters involved in fatigue 
testing under random loading. 
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8 . APPENDICES 


8.1 FORTRAN Programs for Generating Non-Gaussian 
Random Time Histories and Determining 
Certain Peak Statistics 

8.1.1 Binomial Distribution 


C 

C RANDOM TIME HISTORY HAVING A BINOMIAL PROBABILITY 

C DENSITY FUNCTION 

C 

PROGRAM PLOT ( INPUT , OUTPUT , TAPE 21 , TA PE5 =1 NP UT , 

1TAPE6=0UTPUT ) 

DIMENSION P ( 1000 ) , LZ ( 10000 ) , JCOUNT ( 1000 ) , F ( 1000 ) 
DIMENSION LL ( 1000 ),NF( 1000), PI ( 1000 ) 

DIMENSION LZMAX($000 ) , LZMIN { $ 000 ) , ZYM ( 3 ) 

DIMENSION CZMAX( 1000 ),CZMIN( 1000) .XSDBOT(IOO) 
DIMENSION AMID (500) , YJ1 (1000) ,YJ( 1000 ) , YLOGJSM( 500 ) 
DIMENSION P COUNT ( 100 ) , RANDOM (100 ) , YBINTV ( 500 ) 
DIMENSION XMY(1),XXM(1),IN(2),XM(2),YM(2),YYM(3) 
DIMENSION XRMS ( 100 ) , XXMEAN ( 100 ) , XSDTOP ( 100 ) 

DIMENSI ON XAM (U) , YAM ( 2 ) , ZAM ( U ) , YLQGISM ( 500 ) 

C 

C INPUT DATA AND CONSTANTS 
C 

DATA XM/11HP0INT COUNT/ 

DATA YM/13HRAND0M NUMBER/ 

DATA IN/7HHAL SRD, 8 HBINQMIAL/ 

DATA XMY / 9HFREQUEN CY / 

DATA TXM/3.0HCLASS MARK/ 

DATA .YM/2UHMAXIMUM PEAK VALUE COUNT/ 

DATA ZYM/ 2 I 4 HMINIMUM PEAK VALUE COUNT/ 

DATA XAMAOHLOG OF CUMULATIVE FREQUENCY FOR MAXIMUMS/ 
DATA YAM/17 HLOWER CLASS LIMIT/ 

DATA ZAMAOHLOG OF CUMULATIVE FREQUENCY FOR MINI MUMS/ 
ITAPE=6LTAPE21 
C 

C INPUT VARIABLES DEFINED 

C 

C N=NUMBER OF EVENTS 

C M=N UMBER OF RANDOM NUMBERS TO BE GENERATED 

C PEE=PROBABILITY OF SUCCESS IN A GIVEN EVENT 

C BOTTOM=LOWER LIMIT OF RANDOM VARIABLE, LZ 
C TOP=UPPER LIMIT OF RANDOM VARIABLE, LZ 
C FREQM1 =MAXIMUM VALUE OF THEORETICAL AND ACTUAL 

C FREQUENCY OF OCCURRENCE 

C FREQM2=LARGEST PEAK VALUE (MINIMUM) 

C FHEQM3=LARGEST PEAK VALUE (MAXIMUM) 

C 

2 READ (5, 5) N,M, PEE, BOTTOM, TOP, FREQM1 , FREQM2 , FREQM3 
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5 F0RMAT(2I10,6F10.6) 

IF(E0F,5) 991,990 
990 CONTINUE 
AN=N 
P(1)=0.0 
Q=1 . -PEE 
F ( 1 ) =Qm-»N 
KIT=N+1 

WRITE (6,8 ) N,PEE 

8 FORMAT ( 1H1 , 3X, 26HBIN0MIAL DISTRIBUTION N=l£,5X,2EF= 
1 ,F10. 6/ /) 

WRITE(6,19) 

19 FORMAT ( 9X, 1HX, 1 OX , 1HF , II4X , lHP , 12X, 3H1-P/) 

LABCDE=0 

M=M+1 

Jl=2 

J=0 

AM=M 

K0UNT=0 

BINOMIAL PROBABILITY AND CUMULATIVE PROBABILITY 

DO 110 1=1, KIT 

AI=I 

AI1=I-1 

F ( 1+1 ) = ( < AN-AI1 ) /AI )*( PEE/Q )*F (I ) 

P(I+1)=P(I)+F(I) 

P1(I)=1.-P(I) 

IF(P(I+1).LT. 0.00000001) 00 TO 120 
00 TO 30 
120 Jl=I+2 
00 TO 110 
30 CONTINUE 

IF (P(I+1). GT. 0.99999999) GO TO I4.O 
00 TO 110 
14.0 K=I 

GO TO 111 

110 CONTINUE 

111 Jl=Jl-2 
K-K-2 

IF( Jl.EQ.O) GO TO 999 

WRITE (6,1) (I,F(I+l),P(I+2) ,P1(I+1),I=J1,K) 

1 FORMAT (110, 3F15.8) 

00 TO 998 
999 CONTINUE 
KK=K+1 

DO 996 1*1, KK 
11=1-1 

996 WRITE (6, 997 ) II,F(I),P(I+1),P1(I) 

997 FORMAT ( 110, 3F15. 8) 

998 CONTINUE 
J1=J1+1 
K=K+1 
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CALCULATION OP THEORETICAL FREQUENCY OP OCCURRENCE 

DO 112 I=J1,K 
NF(I)=F(I)*AM 
IF(NF(I) .EQ.O ) GO TO 112 
K0UNT=K0UNT+1 
LL(K0UNT)=I 
112 CONTINUE 
MIN*LL(1) 

MAX=LL(KOIJNT ) 

WRITE {6, 63) MIN, MAX 
63 FORMAT ( 5H0MIN=I 6 , 6H MAX=I 6 ) 

GENERATION OF RANDOM NUMBERS 

X=3i£2l637721. 

LSQ=0 
LSUM=0 
53 Y=RANF(X) 

X=0.0 

IF(LABCDE.EQ.M) GO TO $0 
DO 51 I=MIN,MAX 

IF(Y.GE.P(I).AND.Y.LT.P(I+1)) GO TO 52 
GO TO 51 

52 LABCDE=LABCDE+1 
LZ(LABCDE)=I-1 
LSUM=LSUM+LZ (LABCDE ) 

LSQ=LSQ+LZ ( LABCDE )»*2 
GO TO 53 
51 CONTINUE 
GO TO 53 
50 CONTINUE 

XMEAN=FLOAT (LSUM) /FLOAT (M) 

RMS=SQRT (XMEAN* (1 . -PEE )+XMEAN**2 ) 

XMEANSQ=FLOAT ( LSQ) /FLOAT (M) 

STD0EV=3QRT ( XMEANSQ-XMEAH**2 ) 

TIME HISTORY PLOT OF FIRST 100 RANDOM NUMBERS 

MIN=MIN-1 
MAX=MAX-1 
DO 700 1=1,100 
XSDTOP (I ) aXMEAN+ST DDBV 
XSDBOT ( I ) “XMEAN-ST DDEV 
XRMS(I )=RMS 
XXMEAN (I )=XMEAN 
POOUNT(I)*I 
700 RANDOM(I)-LZ(I) 

CALL DDIPLT ( 0 , IN , 100 , PC OUST , XRMS , 1 . 0 , 100 . , BOTTOM, TOP , 
12,XM,2,YM,13,ITAPE) 

CALL DDIPLT ( 0 , IN , 100 , PCOTJNT , XXMEAN ,1.0,100., BOTTOM , 
1T0P,2,XM,2,YM,1U.,ITAPE) 
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CALL DDIPLT ( 0 , IN , 100, PCOUNT , XSDT0P,1 . 0, 100 . , BOTTOM, 
1T0P , 2 , XM, 2, YM , II 4 ., ITAPE ) 

CALL DDIPLT (0, IN, 100, PC0UNT,XSDB0T,1. 0,100., BOTTOM, 
1T0P , 2, XM, 2, YM, II 4 ., ITAPE) 

CALL DDIPLT (1, IN, 100, PCOUNT, RANDOM,!. 0,100., BOTTOM, 
IT OP , 2 , XM, 2 , YM , ll|., ITAPE ) 

KCOUNT-O 

THEORETICAL HISTOGRAM 
WRITE (6,51|.) 

$k FORMAT ( 22H0THE0RETICAL HISTOGRAM//) 

WRITE (6,55) 

55 FORMAT ( 1H0 , 8 X, 1HX, ?X, 11H0CCURRENCES/) 

IF(MIN.EQ.O) GO TO 899 

WRITE ( 6 , 9) (I,NF(I+1),I=MIN,MAX) 

9 FORMAT (110, 115) 

GO TO 898 
899 CONTINUE 
MAX=MAX+1 
DO 896 1=1, MAX 
11=1-1 

896 WRITE(6,897) II,NF(I) 

897 FORMAT (110,115) 

898 CONTINUE 

DETERMINATION OF ACTUAL HISTOGRAM 
WRITE (6,13) 

13 FORMAT (1H0//,3X,16HACTUAL HISTOGRAM//) 

MIN=MIN+1 
MAX=MAX+1 
DO 11 I=MIN,MAX 
IC0UNT=0 
DO 10 J=1,M 

IF(LZ( J) .EQ.I-1 ) IC0UNT=IC0UNT+1 

10 CONTINUE 

11 JCOUNT (I )=ICOUNT 
WRITE (6,55) 

DO 15 I=4GN,MAX 
L=I-1 

15 WRITE (6, 12)L, JCOUNT (I) 

12 FORMAT (110, 115) 

WRITE (6,59) RMS,XMEAN,STDEEV 
59 FORMAT ( 1H1 , 3X, liHRMS^FlO . 3, 3X, 5HMEA1=F10 . 3 , 3X, 15HSTD. 
1DEVIATI0N=F10 . 3 ) 

WRITE (6,56) 

56 FORMAT ( lHO// , 8X , 1HX, 7 X, 1 1HTHE0 . FREQ. ,1|X,12HACTUAL 
1FHEQ. ) 

CHI=0 

LA=MAX 

DO 16 I=MIN,MAX 
L=I-1 



CHI = ( ( (F(I)*AM)-FLOAT(JCOUNT(I)) )**2/(F(I )*AM) )+CHI 

KC0UNT=KC0UNT+1 

AMID(KCOUNT)=L 

YJ1 ( KCOUNT ) =NF ( I ) 

YJ(KCOUNT) =JCOUNT ( I ) 

16 WRITE (6,17) L,NF(I), JCOUNT(I) 

17 FORMAT (110, 2115) 

C 

C PLOTS OF COMPUTED AND THEORETICAL HISTOGRAMS 
C 

CALL DDIPLT ( 0 , IN , KCOUNT , AMI D,YJ1 , BOTTOM, T OP , 0 , FREQM1 , 
11 , XXM, 1 , XMY, 11*, ITAPE ) 

CALL DDIPLT ( 1 , IN , KCOUNT , AMID, YJ,B0TT0M, TOP , 0 , FREQM1 , 1 
1 , XXM, 1 , XMY, 12 , ITAPE ) 

C 

C DETERMINATION OF PEAK AND CUMULATIVE PEAK DISTRIBU- 
C TIONS FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 
C 

MAX=0 

MIN=0 

J=1 

1=0 

200 CONTINUE 
1=1+1 

IF(I.GT.M-l) GO TO 205 
IF(I.NE.l) GO TO 1036 

201 J=J+1 

IF(LZ(I).NE.LZ(J)) GO TO 202 
GO TO 201 

202 CONTINUE 

IF(LZ(I).GT,LZ(J)) GO TO 1037 
GO TO 1038 

1037 MAX=MAX+1 
LZMAX(MAX)=LZ(I ) 

GO TO 1036 

1038 MIN=MIN+1 
LZMIN (MIN ) =LZ (I ) 

1036 CONTINUE 
2OI4. J=J+1 

IF(LZ(I+1).NE.LZ(J)) GO TO 203 
GO TO 201*. 

203 IF(LZ(I+1) .GT.LZ( J).AND.LZ(I+1).GT.LZ(I) ) GO TO 1021 
IF(LZ(I+1) .LT«LZ( J) ,AND»LZ(I+1) ,LT.LZ(I ) ) GO TO 1022 
GO TO 206 

1021 MAX=MAX+1 
LZMAX(MAX)=LZ(I+1) 

GO TO 206 

1022 MIN=4!IN+1 

LZMIN (MIN )=LZ(I+1) 

206 I=J-2 

GO TO 200 
205 CONTINUE 

WRITE (6, 1031*) 
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103k FORMAT (1H1,22X,16HPEAK VALUE COUNT) 

WRITE(6,1039) MAX, MIN 

1039 FORMAT ( 1H0 , 10X, 1I|HMAXIMUM PEAKS=I5 , 2X, lkHMINIMUM PEAK 
lS=l5///kX, 8HINTERVAL, 3X, 3HMAX, 3X, 3KMIN, kX, kHCMAX, i|X, k 
2HCMIN,2X,8HL0G CMAX,3X,8HL0G CMIN/) 

LUINTV =0 
ISUM=0 
JSUM=0 
KC0UNT=0 

1030 LUI NTV =LUI NTV-i-1 

IF ( LUINTV. GT. LA) GO TO 1035 

IC0UNT=0 

JK0UNT=0 

LUINTV=LUINTV-1 

LBINTV=LUINTV-1 

DO 1031 1=1, MAX 

IF(LZMAX(I) .LE.LUINTV. AND.LZMAX(I ) .GT.LBINTV ) ICOUNT= 
1IC0UNT+1 

1031 CONTINUE 
IF(ISUM.EQ.O) GO TO 1000 
ISOM=MAX-ISUM 

GO TO 1051 

1050 ISOM=MAX 

1051 ISUM=ISUM+I COUNT 
IF(ISOM.EQ.O ) GO TO 500 
XL0GISM=AL0G10 (FLOAT (ISOM) ) 

GO TO 510 

500 XL0GISM=0 

510 DO 1032 J=1,MIN 

IF (LZMIN ( J ) . LE. LUINTV . AND. LZMIN (J).GT. LBINTV ) JKOUNT= 
1 JKOUNT+1 

1032 CONTINUE 

IF( JSUM.EQ.O) GO TO 1052 
JSOM=MIN- JSUM 
GO TO 1053 

1052 JSOM=MIN 

1053 JSUM=JSUM+JKOUNT 
IF(JSOM.EQ.O) GO TO 501 
XLOG JSM=AL0G10 ( FLO AT ( JSOM ) ) 

GO TO 511 

501 XL0GJSM=0 

511 WRITE (6, 1033) LBINTV , LUINTV, I COUNT , JKOUNT , ISOM, JSOM, X 
1L0GISM, XLOG JSM 

1033 FORM AT(4l6, 218, F9.k-fFll.k) 

PLOTS OF PEAK AND LOG CUMULATIVE PEAK DISTRIBUTIONS 
FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 

KC0UNT=KC0UNT+1 
YLOGISM(KCOUNT ) =XLOGISM 
YLOGJSM ( KCOUNT ) =XLOG JSM 
YBINTV (KCOUNT ) =LBINTV 
CZMAX(KCOUNT ) -I COUNT 
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CZMIN ( KCOUNT ) =JKOUNT 
AMI D { KC OUNT ) =LUINTV 
LUINTV =LUINTV+1 
GO TO 1030 
1035 CONTINUE 

CALL DDIPLT ( 1 , IN, KCOUNT , YB INTV, YLOGISM, BOTTOM, TOP, 0,1l 
1. ,2,YAM,lt,XAM,12,ITAPE) 

CALL DDIPLT ( 1 , IN , KCOUNT , YB INTV , YLOGJSM .BOTTOM, TOP, 0, it 
1 . , 2 , YAM , it, Z AM , 1 2, IT APE ) 

C ALL DDIPLT ( 1 , IN , KCOUNT , AMID, CZM AX, BOTTOM,! OP , 0 , FREQK 

13.1, XXM , 3 , YYM , 12 , IT APE ) 

CALL DDIPLT { 1 , IN , KCOUNT , AMID, CZMIN ,B0TT0M,T0P , 0 , FRE<31 

12.1, XXM,3,ZYM,12,ITAPE) 

GO TO 2 

99 I STOP 
END 


8.1.2 Poisson Distribution 


RANDOM TIME HISTORY FOR A POISSON PROBABILITY DENSITY 
FUNCTION 

PROGRAM PLOT (INPUT, OUTPUT, TAPE21,TAPE5=INPUT, 
1TAPE6*0UTPUT) 

DIMENSION F( 300 ) ,P ( 300 ) , JCOUNT ( 300) ,LZ (10000 ),LL( 300 ) 
DIMENSION LZMAX ( 5000) , LZMIN ( 5000 ) ,NF ( 300 ) 

DIMENSION CZMAX (1000) , CZMIN (1000 ) ,2YM(3) ,XSDBOT (100) 
DIMENSION AMID(500) ,YJ1 (1000), YJ (1000) ,YL0GJSM(500) 
DIMENSION P COUNT (100) ,RAND0M(100) ,YBINTV (500) 
DIMENSION XMY(1),XXM(1),IN(2),XM(2),YM(2),YYM(3) 
DIMENSION XRMS (100) ,XXMEAN (100) ,XSDTOP (100) 

DIMENSION XAM(lt) ,YAM(2) , Z AM ( it ) , YLOGISM ( 500 ) 

INPUT DATA AND CONSTANTS 

DATA XM/11HP0INT COUNT/ 

DATA YM/13HRAND0M NUMBER/ 

DATA IN/7HHAL 3RD, 7HP0ISS0N/ 

DATA XMY/9HFREQUENCY/ 

DATA XXM/10HCLASS MARK/ 

DATA YYM/2l|HMAXIMDM PEAK VALUE COUNT/ 

DATA ZYM/2itHMINIMUM PEAK VALUE COUNT/ 

DATA XAM/ltOHLOG OF .CUMULATIVE FREQUENCY FOR MAXIMUMS/ 
DATA YAM/17HL0WER CLASS LIMIT/ 

DATA ZAM/ltOHLOG OF CUMULATIVE 7 FREQUENCY FOR MINIMUMS/ 
ITAPE=6LTAPE21 

INPUT VARIABLES DEFINED 

M=NUMBER OF RANDOM NUMBERS TO BE GENERATED 
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U=MEAN OF RANDOM VARIABLE 

BOTTOM=LOWER LIMIT OF RANDOM VARIABLE, LZ 

TOP=UPPER LIMIT OF RANDOM VARIABLE, LZ 

FREQM1 =MAXIMUM VALUE OF THEORETICAL AND ACTUAL 

FREQUENCY OF OCCURRENCE 

FREQM2=LARGEST PEAK VALUE (MINIMUM) 

FREQM3=LARGEST PEAK VALUE (MAXIMUM) 

2 READ (5, 3) M, U, BOTTOM, TOP , FREQM1 , FREQM2 , FREQM3 

3 FORMAT ( 110, 6F10. 6) 

IF(E0F,5) 991,990 

990 CONTINUE 
LABCDE=0 
M=M+1 
N=U99 
Jl=2 
J=0 
AM=M 
K0UNT=0 

IF(U.LE.21. ) XLA-50. 

IF(U.GT.21..AND.U.LE.5l4..) XLA=100. 
IF(U.GT.5iu.AND.U.LE.92.) XLA=l50 . 

IF(U.GT . 92. .AND.U.LE.130. ) XLA=200. 
IF(U,GT.130..AND.U.LE.170.) XLA«2^0. 

IF(U,GT.170. .AND.U.LE.205. ) XLA=300. 

WRITE (6, 102) U,XLA 

102 FORMAT (1H1,3X,29HP0ISS0N DISTRIBUTION, MEAN=F12.7, 
11 OH XUPPER=F6 . 0// ) 

F(1)=EXP (-U) 

P(1)*0.0 

WRITE(6,19) 

19 FORMAT (9X,1HX,9X,UHPR0B,9X,8HCUM PROB/) 

POISSON PROBABILITY AND CUMULATIVE PROBABILITY 

DO 110 1=4, N 
F(I+1)=(U/(I))*F(I) 

P(I+1)=P (I)+F(I ) 

IF(P(I+1).LT. 0.00000001) GO TO 120 
GO TO 30 
120 Jl=I+2 

GO TO 110 
30 CONTINUE 

IF (P(I+1).GT. 0.99999999) SO TO 1*0 
GO TO 110 
kO K=I 

GO TO 111 

110 CONTINUE 

111 Jl=Jl-2 
K=K —2 

IF(Jl.EQ.O) GO TO 999 

WRITE (6,1) (I,F(I+1),P(I+2),I=J1,K) 
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1 FORMAT (110, 2F15. 8 ) 

GO TO 998 
999 CONTINUE 
KK=K+1 

DO 996 1=1, KK 
11 = 1-1 

996 WRITE( 6 , 997 ) II,F(I),P(I+1) 

997 FORMAT (110, 2F15. 8 ) 

998 CONTINUE 
J1=J1+1 
K=K+1 


CALCULATION OF THEORETICAL FREQUENCY OF OCCURRENCE 

DO 112 I=J1,K 
NF(I)=F(I)*AM 
IF(NF(I).EQ.O) GO TO 112 
K0UNT=K0UNT+1 
LL ( KOUNT ) =1 
112 CONTINUE '' 

MIN=LL(1 ) 

MAX=LL( KOUNT) 

WRITE (6, 63) MIN, MAX 
63 FORMAT ( 5H0MIN=I6 , 6H MAX=I6 ) 


GENERATION OF RANDOM NUMBERS 


X=3U521637721. 

LSQ=0 
LSUM=0 
53 Y=RANF(X) 

X=0.0 

IF(LABCDE.EQ.M) GO TO $0 
DO 51 I=MIN,MAX 

IF(Y,GE,P (I) .AND.Y.LT.P(I+1) ) GO TO 52 
GO TO 51 

52 LAB C DE=LABC DE+1 
LZ(LABCDE)=I-1 
LSUM=LSUM+LZ (LABCDE) 

LSQ=LSQ+LZ (LABCDE)«*2 
GO TO 53 
51 CONTINUE 
GO TO 53 
50 CONTINUE 

XMEAN=FLOAT (LSUM) /FLOAT (M) 

RMS=SQRT ( XMEAN* ( XMEAN+ 1 ) ) 

XME AN S Q=FL OAT ( LSQ) /FLOAT (M } 

STDDEV=f3QRT ( XMEANSQ-XMEAN**2 } 

TIME HISTORY PLOT OF FIRST 100 RANDOM NUMBERS 

MIN=MIN-1 
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MAX=MAX-1 
DO 700 1=1,100 
XSDTOP (I ) =XMEAN+STDDEV 
XSDBOT ( I ) =XME AN - ST DDEV 
XRMS ( I ) =RMS 
XXMEAN (I ) =XMEAN 
PCOUNT(I)=I 
700 RANDOM (I ) =LZ ( I ) 

CALL DDIPLT ( 0 , IN , 100, PC OUHT , XRMS , 1 . 0 , 100 . , BOTTOM , TOP , 
12,XM,2,YM,13»ITAPE) 

CALL DDIPLT (0,IN,100,PC0UNT, XXMEAN, 1.0, 100., BOTTOM, 
lT0P,2,XM,2,XM,lli.,ITAPE) 

CALL DDIPLT{0, IN, 100, PC OUHT. XSDTOP, 1. 0,100., ! BOTTOM, 

IT OP , 2 , XM , 2 , YM , 111., ITAPE) 

CALL DDIPLT(0, IN, 100, PCOUHT, XSDBOT, 1.0,100., BOTTOM, 
1T0P , 2 , XM , 2 , YM , 111., ITAPE) 

CALL DDIPLT (1, IN, 100, PCOONT, RANDOM, 1.0, 100., BOTTOM, ' 
IT OP , 2 , XM , 2 , YM , 111., ITAPE ) 

THEORETICAL HISTOGRAM 

KC0UNT=0 

WRITE(6,5lj.) 

$k FORMAT ( 22H0THE0RETICAL HISTOGRAM//) 

WRITE(6,55) 

55 FORMAT ( 1H0 , 8X , 1HX , 7X , 11HOCCORRENCES/ ) 

IF(MIN.EQ.O) GO TO 899 
WRITE(6,9)(I,ffF(I+l),I=MIH,MAX) 

9 FORMAT (110,115) 

GO TO 898 
899 CONTINUE 
MAX=MAX+1 
DO 896 1=1, MAX 
11=1-1 

896 WRITE ( 6 , 897) II,NF(I) 

897 FORMAT (110,115) 

890 CONTINUE 

DETERMINATION OF ACTUAL HISTOGRAM 
WRITE(6,13) 

13 FORMAT ( 1H0// , 3X , 1 6HACTUAL HISTOGRAM//) 

MIN=MIN+1 
MAX=MAX+1 
DO 11 I=MIN,MAX 
IC0UNT=0 
DO 10 J=1,M 

IF(LZ ( J) .EQ.I-1) ICOUNT-ICOUNT+l 

10 CONTINUE 

11 JCOUNT ( I ) = ICOUNT 

WRITE (6, 55) 

DO 15 I=MIN,MAX 
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L-I-l 

15 WRITE(6,12)L,JC0UNT(I) 

12 FORMAT (110, 115) 

WRITE (6, 59) RMS , XME AN , ST DDEV 
59 FORMAT(1H1,3X,4HRMS=F10.3,3X,5HMEAN=F10.3,3X,15HSTD. 

1DEVI ATION=F10 ,3 ) 

WRITE(6,56) 

56 FORMAT (1H0//,8X,1HX,7X,11HTHE0. FREQ * , UX, 1 2HACTUAL 
1FREQ. ) 

CHI=C 

DO 16 I=MIN,MAX 
L=I-1 

CHI = ( ( (F ( I )*AM ) -FLOAT ( JCOUNT (I)) )**2/ ( F ( I ) *AM ) ) +CHI 
KC0UNT=KC0UNT+1 
AMI D ( KCOUNT ) =L 
Y J1 ( KCOUNT ) =NF ( I ) 

YJ (KCOUNT ) =JCOUNT (I ) 

16 WRITE (6, 17 ) L , NF ( I ) , JCOUNT (I ) 

17 FORMAT (110, 2115) 

C 

C : PLOTS OF COMPUTED AND THEORETICAL HISTOGRAMS 
C 

CALL DDIPLT ( 0, IN, KCOUNT, AMID, YJ1, BOTTOM, TOP , 0 , FREQM1 , 
ll,XXM,l,XMY,l4, ITAPE) 

CALL DDIPLT ( 1, IN, KCOUNT, AMID, YJ,B0TT0M,T0P, 0,FREQM1,1 
1 , XXM, 1,XMY,12, ITAPE ) 

C 

C DETERMINATION OF PEAK AND CUMULATIVE PEAK DISTRIBU- 
C TIONS FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 
C 

MAX=0 

MIN=0 

J=1 

1=0 

200 CONTINUE 
1=1+1 

IF(I.GT.M-l) GO TO 205 
IF(I.NE.l) GO TO 1036 

201 J=J+1 

IF(LZ(I).NE.LZ(J)) GO TO 202 
GO TO 201 

202 CONTINUE 

IF(LZ(I).GT.LZ(J)) GO TO 1037 
GO TO 1038 

1037 MAX=MAX+1 
LZMAX (MAX ) =LZ ( I ) 

GO TO 1036 

1038 MIN=MIN+1 
LZMIN (MIN ) =LZ ( I ) 

1036 CONTINUE 
201*. J=J+1 

IF(LZ(I+1) ,NE,LZ( J) ) GOTO 203 
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GO TO 20U 

203 IF(L2 (1+1 ) .GT .LZ ( J) .AND.LZ (1+1 ) .GT.LZ (I ) ) GO TO 1021 
IF (LZ ( 1+1 ) . LT . LZ ( J ) . AND.LZ ( 1+1 ) . LT.LZ (I ) ) GO TO 1022 
GO TO 206 

1021 MAX=MAX+1 

LZMAX (MAX ) =LZ ( 1+1 ) 

GO TO 206 

1022 MIN=MIN+1 

LZMIN (MIN )=LZ (1+1 ) 

206 I=J-2 

GO TO 200 
205 CONTINUE 

WRITE (6,1034) 

1*034 FORMAT ( 1H1 , 22X, 16HPEAX VALUE COUNT) 

WRITE ( 6 , 1039) MAX, MIN 

1039 FORMAT ( 1H0 , 10X, l^HMAXIMUK PEAKS=I5,2X,1UHMINIMUM 

1PEAKS=I5///4X,8HINTERVAL,3X, 3HMAX, 3X , 3HMIN , Ux , UHCMAX, 
24X»4HCMIN,2X,8HL0G CMAX,3X*8HL0G CMIN/) 

LUINT V =0 

ISUM=0 

JSUM=0 

KC0UNT=0 

LA=XLA 

1030 LUINTV =LUINTV+1 
IC0UNT=0 
JK0UNT=0 
LUINTV =LUINTV-1 

* LBINTV =LUINTV - 1 
DO 1031 1=1, MAX 

IF ( LZMAX (I).LE. LUINTV , AND . LZMAX (I) ,GT. LBINTV ) I COUNT= 
1IC0UNT+1 1 

1031 CONTINUE 
IF(ISUM.EQ.O) GO TO 1050 
ISOM=MAX-ISUM 

GO TO 1051 

1050 ISOM=MAX 

1051 ISUM=ISUM+ICOUNT 
IF(ISOM.EQ.O) GO TO 500 . 

XLOGI SM =AL0G1 0 (FLOAT (ISOM) ) 

GO TO 510 

500 XLOGI SM=0 
510 DO 1032 J=1,MIN 

IF ( LZMIN ( J ) . LE . LUINTV . AND . LZMI N ( J ). GT , LBINTV ) JKOUNT= 
1JK0UNT+1 

1032 CONTINUE 

IF( JSUM.EQ.O ) GO TO 1052 

JSOM=MIN-JSUM 

GO TO 1053 

1052 JSOM=MIN 

1053 JSUM=JSUM+JKOUNT 

IF( JSOM.EQ.O ) GO TO 501 
XL0GJSM=AL0G10( FLOAT ( JSOM) ) 
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GO TO 511 
501 XL0GJSM=0 

IF(ISOM.EQ.O. AND. JSOM.EQ.O) GO TO 1035 
511 WRITE (6, 1033) LBINTV , LUINTV , ICOUNT , JKOUNT , I SOM , JSOM , X 
1L0GISM , XLOG JSM 

1033 FORMAT (Ul 6, 218, F9.U.P11.U) 

C » 

C PLOTS OF PEAK AND LOG CUMULATIVE PEAK DISTRIBUTIONS 
C FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 
C 

KC0UNT=KC0UNT+1 
YLOGISM ( KCOUNT ) =XLOGISM 
YLOGJSM ( KCOUNT ) =XLOG JSM 
YBINTV ( KCOUNT ) =LBINTV 
CZMAX ( KCOUNT ) =1 COUNT 
CZMIN ( KC OUNT ) = JKOUNT 
AMID ( KC OUNT ) =LUINTV 
LUINTV =LUINTV+1 
GO TO 1030 
1035 CONTINUE 

CALL DDIPLT (1, IN, KCOUNT, YBINTV, YLOGISM, BOTTOM, TOP, 0,4 

1 . . 2 , YAM , k, XAM, 1 2 , IT APE ) 

CALL DDIPLT(1, IN, KCOUNT, YBINTV, YLOGJSM, BOTTOM, TOP, 0,4 

1.. 2,YAM,4,ZAM,12,ITAPE) 

CALL DDIPLT ( 1 , IN , KC OUNT , AMI D, CZMAX , BOTTOM , TOP , 0 , FRIQM 

13.1, XXM,3,YYM,12,ITAPE) 

CALL DDIPLT ( 1 , IN , KCOUNT , AMI D, CZMIN , BOTTOM, TOP , 0 , FRIQH 

12. 1, XXM,3»ZYW,12,ITAPE) 

GO TO 2 

991 STOP 
END 


8.1.3 Welbgll Distribution 


C 

C RANDOM TIME HISTORY HAVING A WEIBULL PROBABILITY 

C DENSITY FUNCTION 

C 

PROGRAM PLOT ( INPUT , OUTPUT , TAP E21 ,TAPE5=INPUT , 

1 TAPE6=0UTF UT ) 

DIMEN5I ON Z ( 10000 ) , ZMIN ( 5000 ) , Z MAX (5000) , XSDBOT ( 100 ) 
DIMENSI ON XMID ( 100 ) , CZMAX ( 1000 ) , CZMIN ( 1000 ) , ZYM ( 3 ) 
DIMENSION AMID ( 100 ) , YJ1 ( 1000 ) , YJ { 1000 ) , YLOGJSM ( 500 ) 
DIMENSION PC OUNT ( 100 ) , RANDOM ( 100 ) , YBINTV ( 500 ) 
DIMENSION XMY(1),XXM(1),IN.(2),XM(2),YM(2),YYM(3) 
DIMENSION XRMS (100) ,XXMEAN (100 ) , XSDTOP (100 ) 

DIMENSION XAM (Li.) ,YAM( 2) ,ZAM ( U) , YLOGISM (500 ) 

C 

C INPUT DATA AND CONSTANTS 
C 

DATA XM/11HP0INT COUNT/ 
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DATA YM/13HRAND0M NUMBER/ 

DATA IN/7HHAL SRD,10HWEIB,DIST./ 

DATA XMY/9HPREQUENCY/ 

DATA XXM/10HCLASS MARK/ 

DATA YYM/2UHMAXIMUM PEAK VALUE COUNT/ 

DATA ZYM/2UHMINIMUM PEAK VALUE COUNT/ 

DATA XAM/liOHLOG OP CUMULATIVE FREQUENCY FOR MAXIMUMS/ 
DATA YAM/17HLOWER CLASS LIMIT/ 

DATA ZAM/I+OHLOG OP CUMULATIVE FREQUENCY FOR MINIMUMS/ 
ITAPE=6LTAPE21 

INPUT VARIABLES DEFINED 

XINT=CLASS INTERVAL 

ALPHA =PARAMETER OP WEIBULL DISTRIBUTION 
BETA=PARAMETER OP WEIBULL DISTRIBUTION 
ULI MIT =MAXIMUM VALUE OF RANDOM VARIABLE, Z 
M=NUMBER OP RANDOM NUMBERS TO BE GENERATED 
FREQM1 =LARGEST PEAK VALUE (MINIMUM) 

PREQM2=LARGEST PEAK VALUE (MAXIMUM) 

600 READ (5, 601) XINT, ALPHA, BETA, ULIMXT,M,FREQM1 , FREQM2 

601 P0RMAT(UP10.3,I10,2P10.3) 

IP (EOF, 5) 602,603 

603 CONTINUE 
X2=0 
SUMN=0 
M=M+1 

X=31^21637721. 

CHI=0 

AM=M 

WRITE(6,1) ALPHA, BETA 

1 FORMAT (31H1WEIBULL DISTRIBUTION; ALPHA=F10„6,8H 
1BET A=P1 0,6///) 

GENERATION OP RANDOM NUMBERS 

DO 2 1=1, M 
Y=RANF(X) 

■ x=o.o : t*«; : > J 

RWEIB=XWEIB ( Y, ALPHA, BETA ) 

Z ( I ) =RWEIB 
SUMN=SUMN+Z(I) 

2 X2=X2+Z(I)**2 
XMEAN=SUMN/AM 
RMS=SQRT (X2/AM) 

XMEANSQ=X2/AM 

STDDEV=SQRT ( XMEAHSQ-XMEAN**2 ) 

WRITE (6, 1+00) RMS,XMEAN , STDDEV 
Lj.00 P0RMAT(lHl,3X,l4HHMS=F8.l(.,UX,5HMEAN=F8.!4.,3X,l5HSTD. DE 
1VIATI0N=F8.1 j) 
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w 

C TIME HISTORY PLOT OP FIRST 100 RANDOM NUMBERS 
C 

BLIMIT=0 

TOP=ULIMIT 

BOTTOM =BLIMIT 

DO 700 1=1,100 

XSDTOP (I ) =XMEAN +STDDEV 

XSDBOT ( I ) =XMEAN -STDDEV 

XRMS (I )=RMS 

XXMEAN ( I ) =XMEAN i{ 

PCOUNT (I )=I 
700 RANDOM(I)=Z(I) 

CALL DDIPLT { 0 , IN , 100 , PCOUNT, XRMS , 1 . 0 , 1 00 . , BOTTOM, TOP, 
12,XM,2,YM,13, ITAPE) 

CALL DDIPLT ( 0 , IN , 100 , PCOUNT, XXMEAN, 1.0,100. ,BOTTOM, 

IT OP, 2,XM, 2,YM, II 4 ., IT APE) 

CALL DDIPLT(0,IN,100,PC0UNT,XSDT0P,1. 0,100., BOTTOM, 
1T0P , 2, XM, 2, YM, II 4 ., ITAPE) 

CALL DDIPLT ( 0 , IN , 100 , PCOUNT , XSDBOT , 1 . 0 , 100 . ,BOTTOM, 
1T0P , 2,XM , 2,YM, Ilf., ITAPE) , 

CALL DDIPLT (1, IN, 100, PCOUNT, RANDOM, 1 .0 , 100 .^BOTTOM, 
1T0P,2,XM,2,YM,114.,ITAPE) ;• 

C 

C CALCULATION AND COMPARISON OP COMPUTED AND 

C THEORETICAL HISTOGRAMS 

C 

KC0UNT=0 

J=0 ■« 

K=0 

AMDPT=XINT/2. 

WRITE(6,210) > 

210 FORMAT ( 1 H 0 / / , IjX, 10HCLASS MARK, 3X,11HTHE0. FREQ., 2X,12 
1H ACT UAL FREQ.,7X,l4HF(X)/) 

GO TO 7 

6 BLIMIT=BLIMIT+XINT 
IF(BLIMIT.GT.ULIMIT) GO TO 8 
AMDPT=AMDPT+XINT 

J=0 

7 CONTINUE 
DO If. 1 = 1 , M 

IF(Z(I ) .GT.BLIMIT.AND.Z(I) ,LT. (BLIMIT+XINT ) ) GO TO 5 
GO TO k 
5 J=J+1 
Ij. CONTINUE 

ALPHA1 =ALPHA-1 . 

FX=ALPHA*BETA* ( AMDPT**ALPHA1 )*EXF ( -BETA# ( AMBPT&fl'ALPHA 

l)) 

AJ1=FX*XINT*AM 
IF(AJl.EQ.O) GO TO 8 
AJ=J 
J1=AJ1 
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CHI=( ( A J - A J1 ) #*2 ) / A Jl+CHI 
WRITE(6,11) AMDPT, Jl, J,FX 
11 FORMAT (F12.1i,ljX, 110, Ilf>,Fl6. 6) 

PLOTS OF COMPUTED AND THEORETICAL HISTOGRAMS 

K=K+J 

KC OUNT =KC OUNT+1 
AMID(KCOUNT)=AMDPT 
YJ1 (KCOUNT ) =J1 
YJ( KCOUNT )=J 
GO TO 6 

8 CONTINUE . 5 ; 

GALL DDIPLT(0, IN, KCOUNT, AMID, YJ1, BOTTOM, TOP, 0,FREQM1, 
ll,XXM,l,XMY,lfj.,ITAPE) 

CALL DDIPLT(1, IN, KCOUNT, AMID, YJ,B0TT0M,T0P,0,FREQM1,1 
1 , XXM, 1 , XMY, 1 , ITAPE ) 

WRITE(6,10) 

10 FORMAT (1H ///) 

N=M-K 

WRITE(6,10) 

DETERMINATION OF PEAK AND CUMULATIVE PEAK DISTRIBU- 
TIONS FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 

MAX=0 

MIN=0 

M=M-2 

DO 20 1=1, M 
IF(I.NE.l) GO TO 36 
IF(Z(I).GT.Z<I+1)) GOTO 37 
GO TO 38 

37 MAX=MAX+1 
ZMAX(MAX)=Z(I) 

GO TO 36 

38 MIN=MIN+1 
ZMIN(MIN)=Z(I) 

36 CONTINUE 

IF(Z(I+1) .GT.Z(I) .AND.Z(I+1) ,GT.Z(I+2)) GO TO 21 
IF(Z(I+l).LT.Z(I).AND.Z(I+l).LT.Z(I+2)) GO TO 22 
GO TO 20 

21 MAX=MAX+1 
ZMAX(MAX)=Z(I+1) 

GO TO 20 

22 MIN=tGN+l 
ZMIN(MIN)=Z(I+1) 

20 CONTINUE 
WRITE(6,34) 

3k FORMAT ( 1H1 , 22X , 16HPEAK VALUE COUNT) 

WRITE (6, 39) MAX,MIN 

39 FORMAT ( 1H0 , 10X, ll^MAXIMUM PEAKS=l5,2X,ll t HMINIMUM PEAK 
1S=I 5///UX, 8HINTERVAL, 3X. 3HMAX, 3X, 3HMIN, UX, UHCMAX, U", k 
2HCMIN, 2X, 8HL0G CMAX,3X,8HLOG CMIN/) 
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ISUM=0 

JSUM=0 

UINTV=0 

KC0UNT=0 

30, UINTV =UINTV+XINT 
IC0UNT=0 
JC0UNT=0 

BINTV =UINTV-XINT 
DO 31 1=1, MAX 

f IF( ZMAX(I ) ,LE.UINTV.AND.ZMAX(I ) .GT. BINTV) IC0UNT=IC0U * 

XNT+1 

31 CONTINUE 

’ ‘ IF(ISUM.EQ.O) GO TO $0 

I S 0M=MAX- 1 SUM 
GO TO 5l 

50 IS0M=31AX 

51 ISUM=ISUM+ICOUNT 
IF(ISOM.EQ.O) GO TO 500 
XL0GISM=AL0G10 (FLOAT (ISOM) ) 

GO TO 510 

500 XL0GISM=0 

510 DO 32 J=1,MIN 

I F ( ZMIN ( J ) . LE. UINTV . AND. ZMIN ( J ) .GT .BINTV ) JCOUNT=JCOU 
1NT+1 

32 CONTINUE 

IF(JSUM.EQ.O) GO TO 52 

JSOM=MIN-JSUM 

GO TO 53 

52 JSOM=MIN 

53 JSUM=JSUM+JCOUNT 

IF( JSOM.EQ.O) GO TO 501 
XL0GJSM=AL0G1 0 ( FLOAT ( JSOM) ) 

GO TO 511 

501 XL0GJSM=0 

IF(ISOM.EQ.O.AND. JSCM.EQ.O) GO TO 3£ 

511 WRITE(6,33) BINTV , UINTV, I C OUNT, JCOUNT, ISOM, JSOM,XLOGI 
ISM,XLOGJSM 

33 FORMAT ( F6 . 2 , 2H - , Fl*. 2 ,216 , 218 , 2F11 . Ij.) 

PLOTS OF PEAK AND LOG CUMULATIVE PEAK DISTRIBUTIONS 
FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 

KC OUNT=KC OUNT+1 

YLOGISM (KCOUNT ) =XLOGISM 

YLOGJSM ( KCOUNT ) =XLOGJSM 

YBINTV ( KCOUNT ) s®INTV 

IF (KCOUNT . EQ. 1 ) XMID (KCOUNT ) =XINT/2 . 

IF (KCOUNT. GT.l) XMID(KC0UNT)=XMID(KC0UNT-1)+XINT 
CZMAX(KCOUHT )=ICOUNT 
CZMIN (KCOUNT ) =JCOUNT 
GO TO 30 
35 CONTINUE 
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GALL DDIPLT ( 1 , IN ,KCOUNT ,YBINTV ,YLOGISM, BOTTOM, TOP , 0 , 14 . 
1 . ,2, YAM, I 4 . , XAM , 1 , 1 TAPE ) 

CALL DDIPLT ( 1, IN, KCOtJNT, YBINTV,YLOGJSM, BOTTOM, TOP, 0 , h 
1 . , 2 , YAM , k, ZAM , 1 , IT APE ) 

CALL DDIPLT (1, IN, KCOUNT, XMID,CZMAX,BOTTQM, TOP, 0,FREQM 

12 . 1 , XXM, 3 , YYM , 1 , IT APE) 

CALL DDIPLT (1, IN, KCOtJNT, XMID,CZMIN a B0TT0M,T0P,0,FREQM 

12.1, XXM,3,ZYM,1,ITAPE) 

GO TO 600 ' 

602 STOP t ' :? 

END 


8 . 1 . 14 . Exponential Distribution 


RANDOM TIME HISTORY HAVING AN EXPONENTIAL 
PROBABILITY DENSITY DISTRIBUTION 

PROGRAM PLOT (INPUT , OUTPUT , TAPE21 , T APE5=INPUT , 

1TAPE6=0UTPUT ) 

DIMENSION Z ( 10000 } , ZMIN ( 5000 ) , ZMAX ( 5000 ) , XSDBOT ( 100 ) 
DIMENSI ON XMID( 500 ) , CZMAX(1000 ) , CZMIN ( 1000 ) , ZYM (3 ) 
DIMENSION AMID(500) ,YJ1(1000) ,YJ(1000) ,YL0GJSM(500 ) 
DIMENSION PCOtJNT(lOO) , RANDOM (100) ,YBINTV( 500) 
DIMENSION XMY(l) ,XXM(1) ,IN(2) ,XM(2),YM(2) ,YYM (3) 
DIMENSI ON XRMS (100 ) , XXMEAN (100 ) , XSDTOP ( 100 ) 

DIMENSION XAM(l|.) ,YAM(2),ZAM(l4.) ,YL0GISM(500) 

INPUT DATA AND CONSTANTS 

DATA XM/11HP0INT COUNT/ 

DATA YM/13HRAND0M NUMBER/ 

DATA IN/7HHAL SRD,10HEXP. DIST./ . 

DATA XMY/ 9 HPREQUENGY/ 

DATA XXM/IQHCLASSIMARK/ 

DATA YYM/24HMAXIMUM PEAK VALUE COUNT/ 

DATA ZYM/2UHMINIMTJM PEAK VALUE COUNT/ 

DATA XAM/I 4 .OHLOG OP CUMULATIVE FREQUENCY FOR MAXEMUMS/ 
DATA YAM/17HL0WER CLASS LIMIT/ 

DATA ZAM/I 4 .OHLOG OF CUMULATIVE FREQUENCY FOR MINIMUMS/ 
ITAPE=6LTAPE21 

INPUT VARIABLES DEFINED 

M=NUMBER OF UJOM NUMBERS TO BE GENERATED 
THETA=MEAN ANDOM VARIABLE 

FREQM 1 =MAXT VALUE OF THEORETICAL AND ACTUAL 

1FREQUENCY OCCURRENCE 
FREQM2=LA y ^ST PEAK VALUE (MINIMUM) 

XINT=CLAS« INTERVAL 

ULIMIT=UPPER LIMIT OF RANDOM VARIABLE, Z 
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G FREQM3=LARGEST peak VALUE (MAXIMUM) 

600 READ (5, 601) M , THETA , PREQM1 , FREQM2 , XINT , ULtMlT , PREQM3 . 

601 F0RMAT(I10,F10.3»2F10.0,3P10.3) 

IF (EOF, 5) 602,603 

603 CONTINUE 
X2=0 
SUMN=0 
M=M+1 

X=3^21637721. 

CHI=0 

AM=M 

WRITE(6,1) THETA 

1 FORMAT ( 61H1 EXPONENTIAL DISTRIBUTION, F”(X)|THETA*EXP( 

v --- 1=THETA*X) THETA=F10 . 6/ // ) A h 1l: 

C 

C GENERATION OF RANDOM NUMBERS 
C ^ " : 

DO 2 1=1, M 
Y=RANF(X) 

X=0.0 

REXP=XEXP (Y, THETA) 

Z(I )=REXP 
SUMN=SUMN+Z(I) 

2 X2=X2+Z(I)**2 

~ XMEAN =SUMN /AM - - ^ 

RMS=SQRT(X2/AM) 

XMEANSQ=X2/AM 

STDDEV=SQRT ( XMEANSQ-XMEAN**2 ) 

WRITE ( 6 , 14-00 ) RMS, XMEAN, STDDEV ^ 

liOO FORMAT (lHl,3X,UHRMS=F8.U,UX,5HMEAN=F8.1^3Xtl5HSTD. DE 
1VIATI0N=F8 .1}.) 

C .. 

C TIME HISTOHY PLOT OF FIRST 100 RANDOM NUMBERS 

C H- 

J=0 

K=0 

AMDPT=XINT/2. 

BLIMIT=0 

TOP=ULIMIT 

BOTTOM=BLIMIT 

DO $00 1=1,100 

XSDTOP (I ) =XMEAN+ STDDEV 

XS DBOT (I) =XMEAN - STDDEV 

XRMS (I 

XXMEAN(I )=XMEAN 
PCOUNT(I)=I 
500 RANDOM ( I ) =Z ( I ) 

CALL DDIPLT(0,IN,1Q0,PC0UNT, XRMS, 1. 0,100., B0TT0M,TCP, 
12,XM,2,YM,13,ITAPE) 

CALL DDIPLT(0,IN,100,PC0UNT,XXMEAN,1. 0,100., B0TT0M,T0 
IP , 2, XM, 2, YM, l4,ITAPE) 
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CALL DDIPLT (0, IN, 100, PC0UNT,XSOT0P, 1.0,100. , BOTTOM, TO 
' iP,2,XM,2,YM,lMTAPE) 

CALL DDIPLT ( 0, IN, 100, PCOUNT,XSDBOT, 1.0, 100., BOTTOM, TO 
IP , 2 , XM , 2 , YM , llj., ITAPE ) 

CALL DDIPLT ( 1 , IN , 100 , PCOXJNT , RANDOM, 1 .0,100., BOTT OM , TO | 
1P,2,XM,2,YM,114.,ITAPE) ' * 

C: 

' C ' CALCULATION AND COMPARISON OP COMPUTED AND 
C THEORETICAL HISTOGRAMS 
C 

WRITE (6, 210) 

2i0 FORMAT (1HO//,UX,10HCLASS MARK, 3X,11HTHE0. FREQ. ,2X|12?& 
1HACTUAL PREQ.,7X,l4HP(X)/) 

KC0UNT=0 
: GO TO 7 

6 BLIMIT=SLIMIT+XINT 
IP(BLIMET.GT.ULIMIT) GO TO 8 
AMDPT =AMDPT+XINT 

J=0 

7 CONTINUE 

DO k 1=1, M : ? 

" IF(Z(I).GT.BLIMIT.AND.Z(I) ,LT. (BLIMIT+XIHT) ) GO TO $ % 
GO TO k :4 

5 J=J+1 
U CONTINUE 

FX=THETA*EXP ( - AMDPT*THETA ) 

AJ1=FX*XINT*AM 
IF(AJl.EQ.O) GO TO 8 
AJ=J 

J1=AJ1 . 

: : ' * CHI=( ( AJ-AJl )**2)/AJl+CHI 
‘ WRITE(6,11) AMDPT, Jl,J, EX 

11 FORMAT (F12 .l|., 4X, I10,I15,P16 . 6 ) 

:C f i s s ~ 7 ' - ! 

c PLOTS' OF COMPUTED AND THEORETICAL HISTOGRAMS 

C : - . ... 

k=k+j; ; 

’ KC0UNT=KC0UNT+1 
AMID(KCOUNT)=AMDPT 
YJ1(KC0UNT)=J1 * 

YJ(KCOUNT)=J 
GO TO 6 

* 8 CONTINUE J l 

CALL DDIPLT (0, IN,KCOUNT, AMID,YJ1,BOTTOM,TOP , 0, FREQM1 , 
11 , XXM, 1,XMY, lklTAPE) 

CALL DDIFLT ( 1, IN, KCOUNT, AMI D,YJ , BOTTOM, TOP, 0 , FREQM1 , 1 
1,XXM,1,XMY,1,ITAPE) 

WRITE (6* 10)- 
10 FORMAT (1H ///) 

N=41-K 

WRITE (6, 10) 

c ; ; ; ; . . 
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DETERMINATION OP PEAK AND CUMULATIVE PEAK DISTRI- 
BUTIONS FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 
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: . MAX=0 

- MIN=0 

* M=M-2 

DO 20 1=1 ,M 
IF(I.NE.l) GO TO 36 
IF(Z(I).GT.Z(I+1)) GO TO 37 
Si -A GO TO 38 
W-M' MAX=MAX+1 
i& ' T m, ZMAX (MAX ) *Z(I) 

GO TO 36 - 

i: 138 MIN^CN+1 
W -ZMIH (MIH)«Z(I ) 

36 CONTINUE 

U -M IF(Z(I+l).GT.Z(I).AND.Z(I+l).GT.Z(I+2)) GO TO 21 
,IF(Z(I+l).LT.Z(I).AND.Z(I+l).LT.Z(I+2)) GO TO 22 
! .1 GO TO 20 
\:'21 MAX^IAX+1 
if : ■* I- ZMAX (MAX ) =Z ( I +1 ) 

• r GO TO 20 
f 22 MIN=MIN+1 
; • % ZMIN (MIN ) ~Z (I +1 ) 

20 CONTINUE 

WRITE(6,34) 

3U- FORMAT ( 1H1 , 22X, 16HPEAK VALUE COUNT) 

,t WRITE(6,39) MAX, MIN 

39 FORMAT ( 1H0. 10X, 14HMAXIMUM PEAKS=I5, 2X, lUHMINIMUM PEAK 

f 5 is=i5///liX, 8 hinterval, 3X, 3hmax,3x,3Hmin , ux, uhckax, la , a. 

f2HCMIN,2X,8HL0G CMAX,3X,8HL0G CMIN/) 

& -&• ISUM=0 

• '? JSUM=0 

' KC0UNT=0 
UINTV=0 

30 UINTV*UINTV+XINT 
IC0UNT=0 
JC0UNT*O 

- f v BINTVaUINTV-XINT 
i DO 31 1*1, MAX 

*■' - IF(ZMAX(I).LE.UINTV.AND.ZMAX(I).GT.BINTV) ICOUNT=ICOU 
..,,1NT+1 

* 3^ CONTINUE 

IF(ISUM.EQ.O) GO TO 50 
ISOM=MAX-ISUM 
GO TO 51 
$0 ISOM=MAX 
. 51* ISUM=ISUM+ICOUNT 

IF ( ISOM. EQ. 0 ) GO TO 503 
XL0GISM=AL0G10 (FLOAT (ISOM) ) 

GO TO 510 
503 XL0GISM=O 
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510 DO 32 J=1,MIN 

IP ( ZMIN ( J ) . LE . UINTV . AND. ZMIN ( J ) . GT . BINTV ) JCOUNT-JCOU 
1NT+1 

32 CONTINUE / 

IF(J5UM.;||.0) GO TO 52 
J2CM=MXN^isUM 

GO TO >:£#£.- , 

52 JS0M=Ml».4 ' ' 

53 JSUM=J:SDM|JC0UNT 
IP(JSOlj^iO) GO TO 501 
XL0GJSM=&L0G10 ( FLOAT ( JSOM ) ) 

go to mm : 

501 XLOGJSM=Op 

IF(ISO^M.O.AND.JSOM.EQ.O) GO TO 35 
• 511 WRITE(^3ij BINTV, UINTV, I COUNT, JCOUNT, ISOM, JSOM, XLOGX 

1SM,XLQGJSN 

33 F0RMAi|2fi a,2l6,2l8,?9.4,PH. k) 

c &&& - 

C PLOTS -OF. BEAK AND LOG CUMULATIVE PEAK DISTRIBUTIONS 
C FOR BQfH MAXIMUM AND MINIMUM PEAK VALUES 

C - ,-••• 'slf-vJI -• 

KCOUNT=KC'0UNT+1 

YLOGISM ( KCOUNT ) =XLOGISM 

YLOGJSEC KCOUNT ) =XLOG JSM 

YBINTl^feltUNT ) =BINTV 

IF (KCOUNT. EQ.l) XMID(KCOUNT ) =XINT/2 . 

IF (KCOUNT. GT.l) XMID(KC0UNT)=XMID(KC0UNT-1)+XINT 
CZMAX (KGOtiNT ) =ICOUNT 
C ZMIN (KCOUNT ) =JCOUNT 
GO TO 30 '■> 

35 CONTINUE 

CALL DDl-PiT (1, IN, KCOUNT, YB INTV, YLOGISM, BOTTOM, TOP ,0,1*. 

1 . , 2 , yam;©xam, i, itape) 

CALL DDIPLT (1,IN,KC0UNT,YBINTV,YL0GJSM,B0TT0M,T0P,C, u 
1 . , 2 , YAM, !(., ZAM, 1 , ITAPE) 

CALL _DDIPLT(1, IN, KCOUNT, XMID,CZMAX, BOTTOM, TOP, 0,FREQM 
13 , 1 ,XXMy3 , YiM, 1, ITAPE) 

CALL DDIPLT ( 1 , IN , KCOUNT , XMID, CZMIN , BOTTOM, TOP , 0 , FREQM 
12 , 1 , XXMV3^ZYM, 1 , ITAPE) 

GO TO £00- 
602 STOP 
END 



C RANDCMT^ME HISTORY HAVING a LOG NORMAL PROBABILITY 

C DENSITY DISTRIBUTION 

C 

PROGRAM PLOT ( INPUT , OUTPUT , TAPE21 , T APE5=INPUT , 

1TAPE6=0UTPUT ) 
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DIMENSION Z ( lOOOO ),W (10000 ),ZMAX(5000),ZMIN( 5000) 
DIMENSION CZMAX(1000),CZMIN(1000),ZYM(3),XXX(2) 
DIMENSI ON AMID(500)*YJ1 (1000) ,YJ (1000 } ,BMID( $00 ) 
DIMENSI ON PCOUNT ( lQO),, RANDOM ( 100 ) , YBIHTV ( 500) 

DIMENSI ON XMY(l) ,±&(% ),IN(2),XM(2) ,YM(2),YYM(3) , 
DIMENSI ON XRMS (100l*Xf(MEAN ( 100 ) , XSDTOP (100 ) , 

1XSDB0T (100) 

DIMENSION XAM(li),YWg) ,ZAM(l 4 .) ,YL0GISM(500), 

1 YLOGJSM (500 ) 

C ' , 

C INPUT DATA AND CONSTANTS 

DATA XM/llHPOINT 

DATA YM/13HRAND0M NUMBER/ 

' DATA IN/7HHAL SRD/p»OG- NORMAL/ 

DATA XMY/9HFREQUENGYf^ 5 
DATA XXX/17HL0G OF- GLASS MARK/ 

DATA XXM/10HCLASS MM/ 

DATA YYM/ 2 I 4 HMAXIMUMPEAK VALUE COUNT/ 

DATA ZYM/ 2 I 4 HMINIMUM/PEAK VALUE COUNT/ 

- DATA XAM/UOHLOG OF'MbER OF EXCEEDANCES. . . .MAXI MUMS/ 
DATA YAM/17HL0WER^LMS- LIMIT/ 

DATA ZAM/i|.OHLOG 0F s NUMBER OF EXCEEDANCES. .. .MI NIMUMS/ 
ITAPE=6LTAPE21 .. 

C *KTS 

C INPUT VARIABLES DEFINED 

C &. ,vjv. 

C M=’ 3 ER OF RANDOM OTBERS TO BE GENERATED 

C XML-nEAN OF LOGARITHMS OF RANDOM VARIABLE, Z 

C SIGMAS =VARIANCE 0F ; LOGARITHMS OF RANDOM VARIABLE, Z 

C XINT=CLASS INTERVAL # LOGARITHMS OF RANDOM VARIABLE, 
q yz 

C YINT=CLASS INTERVAL/^?/ PEAKS 

C TOP=UPPER LIMIT OF RANDOM VARIABLE, Z 

C BOTTOM=LOWER LIMIT OF, RANDOM VARIABLE, Z 
C BLIMIT=LOGARITHM OF.-^STTOM 

C ULIMIT=LOGARITHM OF TOP 

C FREQM1=MAXIMUM VALUER THEORETICAL AND ACTUAL FREQUE 

C 1NCY OF OCCURRENCE f 

C ' FREQM2»LARGEST PEAK^pLUE (MINIMUM) 

C FREQM3 =L ARGEST PEAK- VALUE (MAXIMUM) 

C * * V 

iq! READ( 5 ,^ 0 ) m,xmu,si®$s,xint,yint,blimit,ulikit,botto 
im,top,freqmi,freqm 2 ;^eqw 3 

lj.0 F0RMAT(I10,2F10.6,2F5.2,UF10.6/3F10.6) 

IF (EOF, 5) 991,990 . •; 

990 CONTINUE . | 

PI =3.110-59265 

AM=M 1 

M=M+1 

X=3l4-521637721. 

BOTTUM=BLIMIT 
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TIP=ULIMIT 
SIGMA=SQRT (SIGMAS ) 

WRITE (6, 2) XMU, SIGMAS, SIGMA 

2 FORMAT ( 8 9H1L0G NORMAL DISTRIBUTION, F(X)=(l/(SIGMA*S 
1QRT ( 2*PI ) ) )*EXP ( - ( LN (X ) -MEAN )**2/ ( 2*SIGMA**2 ) ) //6H ME 
2AN=F5.2,16H SIGMA SQUARE=F£.2,9H SIGMA=F5.2///) 

GENERATION OF RANDOM NUMBERS 


SUMN=0 

X2=0 

DO 1 1=1, M 
Y=RANF(X) 

X=0.0 

RL0GNL=XL0GNL ( Y, XMU, SIGMAS ) 

Z ( I ) =RLOGNL 
W ( I ) =ALOG ( Z ( I ) ) 

SUMN=SUMN+Z (I ) 

X2=X2+Z(I)*»2 
1 CONTINUE 

XMEAN =SUMN/ AM 
RMS=SQRT (X2/AM) 

XMEANSQ=X2/ AM 

STDDEV=SQRT (XMEANSQ-XMEAN**2 ) 

WRITE (6,14-00) RMS, XMEAN, STDDEV 
I 4 .OO FORMAT ( 1H1 , 3X, !|HRMS=F 8 . I 4 ., l|X, 5HMEAN=F8 s I 4 ., 3X, IJjHSTD. DE 
1VIATI0N=F8.14-) 

TIME HISTORY PLOT OF FIRST 100 RANDOM NUMBERS 

J=0 

K=0 

DO 700 1=1,100 
XSDTOP ( I ) =XMEAN+ST DDE7 
XSDBOT(I) =XMEAN - ST DDEV 
XRMS ( I ) =RMS 
XXMEAN ( I ) =XMEAN 
PCOUNT (I )=I 
700 RANDOM (I )=Z(I ) 

CALL DDIPLT ( 0 , IN, 100, PCOUNT , XRMS ,1.0,100. ,BOTTOM,TOP , 
12,XM,2,YM,13,ITAPE) 

CALL DDIPLT (0, IN, 100, PCOUNT, XXMEAN, 1.0, 100 ., BOTTOM, TO 
IP, 2,XM, 2, YM, II 4 ., ITAPB) 

CALL DDIPLT (0, IN, 100, PCOUNT, XSDTOP, 1.0, 100 . , BOTTOM, TO 
1P,2,XM,2,YM,1U,ITAPE) 

CALL DDIPLT (0,IN,100,PC0UNT,XSDB0T,1. 0,100. ,BOTTOM,TO 
IP, 2,XM, 2,YM, ll|.,ITAPE) 

CALL DDIPLT(1,IN,100,PC0UNT,RAND0M,1.0, 100., BOTTOM, TO 
lP,2,XM,2,YM,lq.,ITAPE) 

CALCULATION AND COMPARISON OF COMPUTED AND 
THEORETICAL HISTOGRAMS 
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WRITE(6,210) 

210 FORMAT ( 1H0// , I 4 JC , 7 HMI D- PT . , 2X , 1 1HACTUAL FREQ, 3X, 2l*HL0G 
1MID-PT. THEO. FREQ/) 

KCOUNTsO 
GO TO 7 

6 BLIMIT=BLIMIT+XINT 
IF(BLIMIT.GT.ULIMIT) GO TO 8 ' 

J=0 

7 CONTINUE 
DO k 1=1, M 

IF(W(I ) .GT.BLIMIT. AND«W(I ) .LT. (BLIMIT+.l) ) GO TO ,5 ! 
GO TO 4 
5 J=J+1 
k CONTINUE 

PTL0G=BLIMIT+XINT/2 . 

PT=EXP ( PTLOG ) 

FX=(EXP (- (PTL0G-XMU)**2/(2.*SIGMAS) ) )/SQRT(2.*PI*SIGfM 
IAS) 

J1=FX*XINT#AM 
WRITE ( 6 , 11) PT,J, PTLOG, J1 
11 FORMAT (F10.U., Ill, Fll]., 2, 113) 

PLOTS OF COMPUTED AND THEORETICAL HISTOGRAMS 

K=K+J 

EC OUNT*=EC OUNT+1 
AMID(KCOUNT)=PT 
BMI D ( EC OUNT ) =PTL0G 
YJ1 ( EC OUST ) *=J1 
YJ(ECOUNT)*J 
GO TO 6 

8 CONTINUE 

- CALL DDIPLT (0, IN, ECOUNT, AMID, YJ1, BOTTOM, TOP, 0, FREQMl} 

11, XXM, 1 , XMY, llj., ITAPE) 

CALL DDIPLT ( 1 , IN , EC OUNT , AMID, YJ , BOTTOM, TOP , 0, FREQMl »1 
1,XXM,1,XMY,1,ITAPE) 

CALL DDIPLT (0, IN, EC0UNT,BMID,YJ1,B0TTUM,TIP,0, FREQMl, 

12, XXX, 1, XMY, ll|., ITAPE) 

CALL DDIPLT (l,IN,EC0UNT,BMID,YJ,B0TTIBf,TIP,O,FREQMl, 2 
1 , XXX, 1 , XMY, 1 , ITAPE ) 

N=M-K 

WRITE(6,9) N 

9 FORMAT (1H1,2X,2HN=I5) 

DETERMINATION OF PEA E AND CUMULATIVE PEAK DISTRIBU- 
TIONS FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 

MAX=0 

MIN*=0 

M=M-2 

DO 20 I=1,M 
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( 

IF(I.NE.l) GO TO 36 

IF ( Z ( I ) ,GT.Z(I+1) ) GO TO 37 

GO TO 38 

37 max=max+i 
ZMAX(MAX)=Z(I) 

GO TO 36 

38 MIN=MIN+1 
ZMIN (MIN )=Z(I ) 

36 CONTINUE c 3 

IF(Z(I+l).GT.Z(I).AND.Z(I+l).GT.Z(I+2)) GO TO 2& 
IF(Z(I+l).LT.Z(I).AND.Z(I+l).LT.Z(I+2)) GO TO 22 
GO TO 20 **■■*'* 

21 MAX=MAX+1 

ZMAX(MAX)=Z(I+1) 

GO TO 20 
■' 22 MIN=MIN+1 

ZMIN (MIN )=Z(I+1) 

20 CONTINUE 
WRITE (6, 3k) 

3k FORMAT (1H1,22X,16HPEAK VALUE COUNT) 

WRITE( 6, 39 ) MAX,MIN« . ; 

39 FORMAT ( 1H0.10X, UjHMAXIMOH PEAKS=I5, 2X, lUHMIHIMUft PEAK 
1S*I5//AX, 8HINTEHVAL, 3X-3HMAX, 3X, 3HMIN, IjX, IjHCMAf , UX, k 
2HCMIN,2X,8HLOO CMAX,3X,8HE0G CMIN/) 

C ** COUNTS PEAK VALUES ** 

UINTV =0 

- • ISUM=0 ' 

JSUM=0 

KC0UNT=0 

30 UINTV =UINTV+YINT 
IC0UNT=0 
JC0UNT=0 

BIKTV=UINTV-YINT 

PTL0G=SINTV+YINT/2. 

DO 31 1=1, MAX ‘ 

IF(ZMAX(I ) .LE.UIHTV. AND.aiAX(I ) .GT.BINTV) IC0UNT=IC0U 
1NT+1 

31 CONTINUE 
IF(ISUM.EQ.O) GO TO 50 
ISOM=MAX-ISUM 

GO TO 51 

50 ISOM^MAX 

51 I SUM=I S UM+I C OUNT . ; . 

IF(ISOM.EQ.O) GO TO 500 
XL0GISM=AL0G10 (FLOAT (ISOM) ) 

go T 0 510 

500 XL0GISM=0 
510 DO 32 J=1,MIN 

IF ( ZMIN ( J ) .LE . UINTV . AND. ZMIN ( J ) .GT.BINTV ) JC0UN$4=JC0U 
1NT+1 

32 CONTINUE 
IF(JSUM.EQ.O) GO TO 52 
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JSOM=MIN- JSUM 
GO TO 53 

52 jsom=min 

53 JSUM=JSUM+JCOUNT 
IP(JSOM.SQ.O) GO TO 501 
XL0GJSM=AL0G10 (FLOAT ( JSOM ) ) 

GO TO 511 

501 XL0GJSM=0 

IF(ISOM.EQ.O.AND.JSOM.EQ.O) GO TO 35 
511 WRITE ( 6 , 33) BINTV,UINTV,ICOUNT,JCOUNT,ISOM,JSOM,XLOGI 
1SM,XL0GJSH 

33 F0RMAT(2F6. 1,216, 218, F9.U,F11.W 
C 

C PLOTS OF PEAK AND LOG CUMULATIVE PEAK DISTRIBUTIONS 
C FOR BOTH MAXIMUM AND MINIMUM PEAK VALUES 

Q 

KG OUNT =KC OUNT+1 
YLOGISM (KCOUNT ) =XLOGISM 
YLOGJSM (KCOUNT ) =XLOG JSM 
YBINTV (KCOUNT ) =43INTV 
C ZMAX (KCOUNT )=IC OUNT 
CZMIN (KCOUNT )=JC OUNT 
AMID (KCOUNT ) =PTLOG 
GO TO 30 
35 CONTINUE 

CALL DDIPLT(1, IN, KCOUNT, YBINTV , YLOGISM, BOTTOMiTOP, 0,14.' 
- 1, ,2,YAM,4»XAM,1,ITAPE) 

CALL DDIPLT ( 1 , IN , KCOUNT , YBINTV , YLOGJSM , BOTT OM , T OP , 0 , 14 . 
l.,2,YAM,U.,ZAM,l,ITAPE) 

CALL DDIPLT(1, IN, KCOUNT, AMID, CZMAX,B0TT0M,TQP,0,FREQM 
13,1,XXM,3,YYM,1,ITAPE) 

CALL DDIPLT (1 , IN, KCOUNT, AMID, CZMIN, BOTTOM, TOP , 0 , FREQM 
12,1",XXM,3,ZYM,1,ITAPE) 

GO TO ui 
991 STOP 
END 



8.2 FORTRAN Program Subroutines 
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8. 2. 1.2 Plotting Subroutine DDIPLT 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 

G 

C 

C 

C 

G 


IN 

N 
X 
Y 
XMIN 
XMAX 
YMIN 


NXM 

XM 

NYM 

YM 


y 


SUBROUTINE DDIPLT (IEC,IN,N,X,Y,XMIN,XMAX,YMIN,YMAX, 

L NXM,XM,NYM,TM,ISYM, DDITPE) 

I EC DISPLAY END CODE 

0 DATA INCOMPLETE 

1 DISPLAY COMPLETE . ^ * 

PROGRAMMER INITIALS AND JOB IDENTIFICATION { 

2 WORD ARRAY) ? , r 

NUMBER OF DATA POINTS TO BE PLOTTED 
ARRAY CONTAINING X COORDINATES 
ARRAY CONTAINING Y COORDINATES 

*** DATA RANGE TO BE PLOTTED. IFVALUES 
*** ARE ZERO THE RANGE WILL BE DETERMINED 
*** ON FIRST CURVE OF DISPLAY 
YMAX ss-sh* 

NUMBER OF WORDS IN HORIZONTAL MESSAGE ARRAY 
HORIZONTAL MESSAGE 'T 

NUMBER OF WORDS IN VERTICAL MESSAGE^ AR&AY 
VERTICAL MESSAGE 
ISYM SYMBOL CODE 

DDITPE TAPE NAME IN PROGRAM CARD V,: r, 

''’I* .. 

DIMENSI ON IN ( 2 ) , ID ( 3 ) , INIT ( 9 ) , PLOT (100), IPLOT ( 100 ) 
DIMENSION X(N ) , Y(N), ISTAB(llj.) a SM(3) 

DATA (ISTAB(I),I=l,li|.)/076203655337575757575 ? 

1 07620365531^757575757# -4 

2 076203655207575757575, 

3 07620365553757575757# , U 

k 07620365554757575757# > 

5 076203655U-07575757575, ; 

6 076203655617575757575* - 

7 076203655157575757575, 

8 076203655167575757575, 

9 076203655137575757575* * I 

1 076203655111-7575757575, 

2 076203655327575757575. 

3 076203655527575757575, 

4 ■ 0766055757575757 57 575/ Y 

data (SM(l),l=l,3)/077757575757575757575,075757575757 
1575757575, 075757575757575757575/ . p , 

DIMENSION XGRID(lli-), YGRID(ll4.),IXGRID(lli),IYGRID(lU) 
DIMENSION XSC (II 4 .) ,IXPOS (ll|.) ,YSC (Ik) ,IXSC(lli) ,IYSC(ll}.) 
EQUIVALENCE (PLOT,IPLOT) 

EQUIVALENCE (XGRID,IXGRID) , (YGRID,IYGRID) 

EQUIVALENCE (XSC , IXSC ) , (YSC , IYSC ) . 

INTEGER DDITPE 

DATA NP AGE/O/, IEND/l/,NCAM/076002275757575757575/ 

DATA (ID(I ) , 1=1 , 3 )/07 600227654123717377 5 , 00,075757575 
17575757600U4/ 
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- DATA (INIT(I ), 1-1,9) /07 6)4>3U1337377 57 57575, 00 , 076U.0^. 

113371775757575.00. 075757575757575757514), 00,076U03Ul33 

261775757575. 00. 075757575757575760014;/ 
c 

C TRANSFER ON FIRST ENTRY 
C 

IF (INIT(2).EQ.O) GO TO 100 
C 

C TRANSFER FOR START OF NEW DISPLAY (IEND=1) OR SET 
C 1IG0DE = 1 
C 

IF (IEND.EQ.l) GO TO 110 
K=1 

GO TO 1U0 
G 

C SEND INITIAL IDENTIFICATION FRAME 

C 

100 CALL DISBCD (IN(1) ,INIT(2) ,1) 

ID(2)=INIT(2) 

WRITE (DDITPE) ID 
G 

C INITIALIZE FOR A DISPLAY 

G 

110 NPAGE=NPAGE*1 
LPCT=0 

IPLOT(l) = NOAM 
K=2 

GALL DISBCD (IN(2) ,INIT(1*),1) 

G 

C CHECK FOR X MAX-MI N VALDES 
C 

IF (XMIN.EQ.0.0.AND.XMAX.EQ.0 o 0) GO TO 120 
XMN=XMIN 
XMX=XMAX 
GO TO 125 
120 XMN=X(1) 

XMX=X(1) 

DO 122 1=2, N 

IF (X(I).LT.XMN) XMN=X(I) 

IF (X(I).GT.XMX) XMX=X(I) 

122 CONTINUE 
C 

, C CHECK FOR Y MAX-MIN VALUES 

U 125 IF (YMIN .EQ. 0 . 0 . AND.YMAX. BQ. 0.0) GO TO 12? 

YMN=YMIN 
- YMX=YMAX 

GO TO 130 
127 YMN=Y(1) 

YMX=Y ( 1 ) 

DO 129 I=2,N 

IF (Y(I).LT.YMN) YMN=Y(I ) 
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c 

c 


c 

c 

c 


c 

c 

c 


•» 

#> 

* 

* 


c 

c 

c 


c 

c 

c 


. IP (Y(I).GT.YMX) YMX=Y ( I ) 

129 CONTINUE 

^ETUP HORIZONTAL (X) MESSAGE 

*i fii' ■ 

1S0: If os=520-nxm*Uo 

*f;.pIGH=IPOS/UOB 
' ° LOW=I P OS -IHIGH#i|.OB 

p^^LOT(K)=7575757576l|-OOOOOOOOOB+(IHIGH*100B+ILOW)*xOGO 

I&#K+1 

t^CALL DISBCD (XM,PLOT(K),NXM) 
l?ftKt%K+NXM 

TUP VERTICAL (Y) MESSAGE 


P0S=520-NYM^0 
-V* t%igh=IPOSAOB 


J^^dW»IPOS-iHIGH*^OB 


imL DISBCD (YM,PLOT(K),NYM) 
^|pK+HYM 


MINIMUM VALUES FOR DISPLAY 


ENCODE (10,1,XMNT) XMN 
d WR HAT (E10.3) 

|f:#CODE (10,1,YMNT) YMN 
a 1 : -M TO 150 

tullDlTIONAL DATA FOR DISPLAY 
wjl&ECK X MAX-MIN VALUES 

ll*.0 IF (XMIN.NE.O.O.OR. XMAX.NE.O.O) GO TO li*2 
• -- IC0DE2=0 

^TO 11+5 

li& ,Xp*XMIN 
XKX=XMAX 
•Vf)S0DE2=l 

CHECK Y MAX-MIN VALUES 

•* • - * 

l)+5 1ft (YMIN.NE.O.O. OR. XK.1X.NE, O.O) GO TO 11*7 
^ $F (IC0DE2) 150,124.7,150 
IU.7 YMN=YMIN 
YKX=YMAX 

^$UECK DATA RANGE, ADJUST IF ZERO 

150 XRANGE=XMX-XMN 
YRANGE=YMX-YMN 
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C 

C 

C 


IP (XRANGE.NE.O.O) GO TO 154 
IP (XMN.E^.d.O) GO TO 152 
XMN=XMN- . ltfABS ( XMN ) 

XMX=XMX+ . 1MBS (XMX) 

GO TO 154 M 

XMN=-1 . 0 
XMX=1.0 ;'r>Sg. 

IP (YRANGE^ffi.O.O) GO TO 160 

ip (ymn.I&Mo) go to 156 

YMN=YMN- M*mS (YMN ) 

YMX=YMX+ .1*|BS(YMX ) 

GO TO l6p>§§, , 

nor— x # og|^L 

YMX=1.0 MM .. 

FIND INCifi®T AND ADJUST MAJ 


IT AND ADJUST MAX-MI N VALUES 


)Ao -° 

IP (XINCafeMo.O) GO TO 164 

1 = 1+1 

XING=XING:pOSO 

go to lemjlt 

IP (XINC Mo) GO TO 166 

1 = 1-1 

XINC=XING-»10;0 
GO TO 164 

IP (XINC Arl.O) GO TO 116? 
IF (XINC.ItJI. 0) GO TO 16? 
XINC=2.0 A . g. 

GO TO I69f, m.. 

XINC=1 . 0 w 
GO TO 169 »>'%*••• 

IP (XINC.GTi5iO) GO TO 168 
xinc=5.o 
GO TO 16% .r-U ♦ 
XINC= 10 . 0 'T^' 

IIX=I 

IP (I) 171 ^ 3,172 

XMN=XMN*lOr 0 ^ 

XMX=XMX*1 0.0 
GO TO 170W - 

1=1-1 

xiNc=xiNc»ro;o 

GO TO 170 - 

XR=AINT ( XMN/2TNC ) *XINC 

IP (XR.GT~XM&) XR=XR-XINC 

XMN=XR 


XR=AINT (XMX/XINC )*XINC 
IP (XMX.GT.XR) XR=XR+XINC 
XMX=XR 
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YINC = (YMX-YMN)/10.0 

1=0 \V 

Ilk IP (YINC.LT. 10.0) GO -TO 175 

1=1+1 r *• 

YINC=YINC/10.0 
GO TO 1714- 5 

175 IF (YINC.GE.1.0) GO TO 1| 6 


k 4 / 




+^3, .,»*•’ •' 


1=1-1 

YINC=YINC*10.0 
GO TO 175 

176 IP (YINC.EQ.1.0) GO 
IP (YINC.GT.2.0) GO 
YINC=2.0 
GO TO 179 

1177 YINC=1 .0 
GO TO 179 

177 IP (YINC.GT.5.0) GO 
YINC=5.0 
GO TO 179 

178 YINC =10.0 

179 HY=I 

180 IP (I) 181,183,182 

181 1 = 1+1 
YMN =YMN*10 . 0 
YMX=YMX*10 . 0 
GO TO 180 

182 1 = 1-1 

YINC =YI NC*10 . 0 
GO TO 180 

183 YR=AINT ( YMN/YINC )*YI&0 1 
IP (YR.GT.YMN) YR=YR-tENf 

YMN=YR : y 

YR=AINT ( YMX/YINC )*Y£HC 
IP (YMX.GT.YR) YR=YR+YINC 
YMX=YR - I*. 

1 - • * y * 

! TABLE FOR X AND Y GHJO .... 


wA ? 

■ 

2%: Mi' 

£> V 

y #>■ 

;• ;<&■; yy ' 

‘"'Q; 

i-S&I 's&k. 

£ a- 




m 


. *-> 


XSCALE=( XMX-XMN )/983l0 
YSCALE=( YHX-YMN )/98340 
IP (IEND.EQ.O) GO TO 21*6 
IX=1 

XGRID(1)=XMN ■ ^ 

I8I4. IX=IX+1 

XGRID(IX)=XGRID(IX-1)+XINC 
IP (XGRID(IX).LT.XMX):GO TO 184 

IY=1 i ; ’ 

YGRID(1)=YMN 
186 IY=IY+1 

YGRI D (I Y) =YGRID(IY-1)+YINC 
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13a 


c 

c 


IF (YGRID(IY).LT.YMX) GO TO 186 
TABLE FOR X AND Y LABELS 


■1 

■ ... t. 

k 

•- ■ 


y 

; 


% 


DO 190 1=1, IX 
190 IXSC (I )=XGRID(I ) 

192 IF (IABS(IXSC (1) ) .LT.iOOO„AND.IABS(iX&C (IX) ) .LT.IOOO) 

1GO TO' 196 ' 

DO 193 1=1, IX 

193 IXSC ( I )=IXSC (I )/10 
GO TO 192 

196 DO 205 1=1, IX 
IF (IABS (IXSC (I ) ) .GT.9) GO TO 197 
IXPOS (I )=-2lj. 

GO TO 200 

197 IP (IABS (IXSC (I)).GT. 99) GO TO 
IXPOS (I )= -20 
GO TO 200 

198 IXPOS (I )=-l 6 
200 ENCODE (10,2,IXSC(I) )IXSC(I ) 

2 FORMAT 
205 CONTINUE 

IXPOS (IX )=-2li 

CALL DISBCD (IXSC(l) ,IXSC(1) ,IX) 

DO 210 1=1, IY 
210 IYSC (I )=YGRID(I ) 

212 IF (IABS(IYSC(1)). LT.IOOO. AND. lABS^mC(lY)). LT.IOOO) 
1G0 TO 216 

DO 213 I=1,IY 

213 IYSC (I )=IYSG (I )/lO 
GO TO 212 

216 DO 220 1=1, IY 

ENCODE (10, 2, IY&C (I ) )IYSC (I ) 

220 CONTINUE 

CALL DISBCD (IYSC(1),IYSC(1),IY) 

do 230 1 = 1 , ix ; yi 

230 IXGRID(I )=(XGRID(I )-XMN )/XSCALE +^ '$ - * 

DO 235 I=1,IY • m ** 

235 IYGRID(I )=(YGRID(I )-YMN)/YSCALE+32.W 


pr, m. 

$* /;; . 

Ki- ' - 

f-r- >!*■*'■ 

dM ■ 


I w 


■fl ' Wv 
*• V-i- • 

;'?f . \ :<*•• 

'ft/' 


DRAW GRID LINES 


T- 


238 iai=ixgrid(i)Aob r/tv 

IA2=IXGRID(1 )-IAl*40B 
IA3=IXGRID(IX)A0B * 

IAl4.=IXGRID(IX)-IA3*l».OB >• * 

IA5=76600000000000000000B+(IAl*100fe+JA2)*100000000000 
10B+ ( I A3*100B+I Al|.)*10000B f 

I A6=7 640001U000075757575B v ! • 

DO 2l}.0 1=1, IY 
IA1=IYGRID(I ) AOB 
I A2=IYGRID(I )-IAl*UOB 
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IA3=IA1*1003+IA2 
IPLOT (K)=IA5+IA3<'100000000B+I A3 
I PLOT (K+l J =IA6+IA3*100000000B 
IPLOT ( K+2 ) =1 YSC { I ) 

2LuO K=K+3 

I A1=I YGRID ( 1 ) /I 4 .OB 
I A2=I YGRID ( 1 )-IAl*40B 
IA3=IYGRID(IY)AOB 

I A4=I YGRID ( I Y) -I A3 AOB §■ , . 

IA5=766OOOOOOOOOOOOOOOOOB+(IA1#1OO3+IA2)‘&1OO0qC>00OB+( 
1IA3*100B+IA4) *" 

IA6=76l4.00000001775757575B 
DO 245 1=1, XX 
IA1=IXGRID(I ) AOB 
IA2=IXGRID(I )-IAl#!+OB 
IA3=IA1*100B+IA2 
IA7=IXGRID(I I+IXP03 (I ) 

IA8=IA7/40B 

IA9=IA7-IA8AOB 

IA10=IA8*100B+IA9 

IPLOT (K)=IA5+IA3*1000000000000B+IA3*10000B 
IPLOT (K+l ) =IA6+IA10»1000000000000B 
IPLOT (K+2)=IXSC (I ) 

245 K=K+3 

2 I 4.6 IP (IIX.GE.O) GO TO 2 I 4.7 
IIX=IIX+1 
’ XMN=XMN/10.0 
XMX-XMX/lO . 0 
XSCALE=XSC ALE/10 . 0 
XINC=XIN C/10,0 
GO TO 246 

247 IP (IIY.GE.O) GO TO 1180 
IIY=IIY+1 
YMH=YMN/10.0 
YMX=YMX/10.0 
YSCALE=YSCALE/lO . 0 
YINC=YINC/10.0 
GO TO 247 

C 

C DISPLAY MINIMDM VALUES AND INCREMENT 

c - 

1180 IP (IEND.NE.l) GO TO 1247 

I PLOT ( X ) =7 640002437 3775756060B 
IPLOT (K+l ) *606067604431456013608 
CALL DISBCD (XMNT, IPLOT (K+2), 1) 

IPLOT (K+3 ) =60314523512544254563B 
1185 ENCODE ( 10,1, IPLOT (K+4)) XINC 

CALL DISBCD (IPLOT (K+4), IPLOT ( K+4), 1) 

IPLOT (K+5)=60607060443145601360B 
CALL DISBCD ( YMNT, IPLOT (K+6),l) 

IPLOT (K+7 ) =603145235125442545638 

1195 ENCODE (10,l,IPL0T(K+8)) YINC 
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CALL DISBCD (IPL0T(K+8),IPL0T(K+8),1) 

K=K+9 

JSET SYMBOL AND SCALE DATA TO PLOT 

1247 1?L0T(K)=I3TAB(ISYM) 

- K=K+1 

: >4X) 250 1=1, N 
• tP(X(I ) .GT.XMX) GO TO 301 
j^P(X(I ) .LT.XMN) GO TO 302 
. ;f P ( Y ( I ) . GT . YMX ) C-0 TO 303 
. |F(Y(I).LT.YMN) GO TO 3 OI 4 . 
j©5 ; 2A1=(X (I ) -XMN )/XSCALE+I;0 .5 
f V; A2=IAl/l4.0B 
iu A3=I Al-I A2*40B 

v?MiAU»(Y(I )-YMN)/YSCALE+32.5 
s ! -■/- |1a5 *1 aU/ ioB 
'IA6=IA4-IA5*40B 

f IPLOT ( K ) =7575757^757^00000000B+ (IA2*100B+I A3 )*10000B+ 
HA5»100B+IA6 
U t=K+l 
'GO TO 2^9 

3Gl'fP(X(I).EQ.O.) GO TO 2l|S 
p||P(ABS((X(I)-XMX)A(I)).OT.5.E-ll4.) GO TO 248 
fO TO 305 

302 tP(X(I).EQ.O.) GO TO 248 
:| t ^|F(ABS((XMN-X(I))/X(I)).GT.5.E-ll|.) GOTO 2U8 
iV @0 TO 305 

303 IF(Y(I).EQ.O.) GO TO 21$ 

. |F(ABS((Y(I)-YMX)/Y(I)).GT.5.E-14) GOTO 248 
M ? '#0 T0 305 

304. fF(Y(I).EQ.O.) GO TO 248 

|f(abs((ymn-y(i))A(i)).ot.5.e-iU) go to 248 
I ; |0 TO 305 
248:iiPCT=LPCT+l 
249 YF (K.LT.101) GO TO 250 

WRITE (DDITPE) (IPLOT( J), J=l,100) 

25d 3 'C0NTINUE 

' TF (K.EQ.l ) GO TO 252 
'£=K-1 

WRITE (DDITPE) (IPLOT(I ),I*1,K) 

252-IP (IEC.EQ.l) GO TO 255 
f END=0 
GO TO 1000 

255 Encode ( 10 , 2 , init ( 6 ) jnpage 
Gall disbcd (init(6),ihit(6),d 
iend=i ■ 

IF (LPCT.EQ.O) GO TO 2$6 
ENCODE (10,2,INIT(8))LPCT 
CALL DISBCD (INIT(8),INIT(8),1) 

GO TO 260 
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256 INIT ( 8 ) =60606060606060606060B 
260 WRITE (DDITPE) (INIT (I), 1=1,9) 
WRITE (DDITPE) (SM(I ), 1=1, 3) 
IEND=1 

1000 RETURN vV? 

END « v*. 

i'.' 'ii: 


8.2.2 Funct Ion ;g^.bpr ograma 

/ i-, 

8,2. 2.1 Fmpltlon XWEIB 


FUNCTION^XfElB (FX, ALPHA, BETA) 
ALPHA1*&#&PHA 
XWEIB =( (M1.%ETA)*AL0G(1.-EX) )*»(ALPHA1) 
RETURN ' "” v ' 

END 



8 . 2 . 2. 2 

w~- 

(FX, THETA) 
XEXP=XWfflf SX, 1 . , THETA ) 
RETURN rffcSk 
END 2$ 'C 

WA, , 

8. 2. 2. 3 function XLOQNiu 


FUN C TI ONv&gGNL ( FX , XMU, SIGMA ) 

SUM1=FX '4mm ' 

IF(SUM1.#.J5) SUM1=1.0-SUM1 
XNU=SQRT$U*lG( 1 . 0/ (SUM1*SDM1 ) ) ) 

Z=XNU- ( (Mi5^5l7+.802853*XNU+.010328*XNU»XNU)/(1.0+1. 
114.32788*XNU#^189269*XNU#XNU+. 001308*XNU«t3 .0) ) 
IF(FX.LTv.5tZ=-Z 
XLOGNL=EXP ( SQRT (SIGMA )*Z+XMU) 

RETURN 

END 
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8.3 Typical FORTRAN Program Block Diagram 
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CF{x) 

FU) 

MEAN 

m 

n 

P 

q 

RMS 

X,Z,LZ 

Y 

a 

P 

r( ) 

9 

p 

^log 

a log 


8.U List of Symbols 

- Cumulative probability density function 

- Probability density function 

- Arithmetic mean of random variable 

- Integer 'Vi 

- Number of events (positive integer) 

- Probability of success 

- Probability of failure 

- Root -mean-square 

- Random variables 

- Uniformly distributed random number 

- Parameter of Weibull probab±i||| density function 

- Parameter of Weibull probability density function 

- Gamma function 

- Parameter of Exponential probability density 
function 

- Parameter of Poisson probability density function 

- Mean of logarithms ' 

- Standard deviation of logarithms 
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8.5 Derivation of Theoretical Equations Describing 
Properties of Density Functions 

8.5.1 General 

Mean: E(x) = j xF(x)dx 

U — OO 

Mean square: E(x 2 ) = j x 2 F(x)dx 

U*» 

Variance: a 2 = E(x 2 ) - j~E(x)J 2 

8.5.2 Exponential Amplitude Probability Density Function 

Probability density function: F(x) = 0 exp(- 

Mean: E(x) = j 

Jo 


& 


x 0 exp(- 0 x)dx = ~ 
0 




Mean square: E(x 2 ) 


2 1 

Variance: a = — 5 


■£ 


x 2 0 exp(- 0 x)dx = 


T. 

tit 






8.5.3 Welbull Amplitude Probability Density Fun^ilojit - 

f*:$V 

Probability density function: F(x) = a0x a "*exp(-0x a ) 

Mean: E(x) = a 0 x a exp { - 0 x a } dx 


Let 


then 


and 


0 x a = z 


( Z Y^ / a 

(?) “ X 


apx 


dz 

all 


dx * 



1U1 


Substituting we get 


E(x) = P az exp (-z) — ~ — 
dO _ Q _a-l 


Simplifying we get 


E(x) 


I^ /a P" 


= (tr £ 


(z) exp(-z)dz 


Prom Ryshik and Gradstein (I 963 ) equation (3.225) we get 


-w - W) 


Mean square: 


E(x 2 ) = a(Jx a+ ^exp(-Px a )dx 


z = px® 


then 


dx = 


Substituting and simplifying we get 


E(x 2 ) 


z*'®exp(-z)dz 


Prom Ryshik and Gradstein (1963) equation (3.225) we get 

» ! > ■ (if 



1U2 


Variance: c 2 = (|) ' - [''(HrO ^ 

• k 
V 

8.5»liy D§g-Normal Amplitude Probability Density Function 


V / £ 

Probability density function; 


^'4 ’ ' ' w 


Let 


JH- v -w . 

x ) = 

pit': • .'V:S^V 

SS'Vj^ •; * 

£. 3 *- X^wr ' 


Xff. /2it 
log 


exp 


2 a 


2 (^°ge x ” ^log) 
r log 


then 


' ** V V*' * * 

rj,:; - 3^; 

>*-»„•' 0 •* 
* v ' ; *4# 

&* 

••■tJvT' JV***', 


z = log e x 


P(x) = G(log e x) a G(z) 

E(x n ) » P x n dP(x) 
f UO 


Substituting we get 

E(x n ) = 


m'm 

TJtfK ’ J-g 






exp(nz)dG(z) 


Prom Ait chi son and Brown (1957) we get 


4>* * *<f. 
< * • 

•* 




E(r n ) = «qp(n|l log + 1 n 2 cf og ) 


Therefore 


Me^ E(x) *, expfitiog + | ®f 0 g) 

Mean square; E(x 2 ) = exp^2y, lQg + 2<y fog) 



Variance: c 2 = «p(2n log + 2c 2 og ) - «p(2|* log * c 2 0g ) 


Simplifying we get 

cijvi' 


02 - Ift^log + °lcg 

• ■‘* 1 % . 

&ty.m 


)R»U - *] 


8 Poisson Amplitude Probability Density Function 

-V 1 X 

Probability|delsity function: F(x) -=~ e - ■ ^ 

•*i**-, -*fY' 

xe ^|x 

. xl 

*?A: P 

But at x = 0 ttfis Jirst term of the series is zero, 

m.$y 

therefore * 


Mean: E ( x } 


JL --m.x ® -U..X-1 


Let z = x - 1, : the& 

:"&■ '-m- 

■> ••:•% •. sf* 


«*• ' V 




oo —UL Z 

E(x) = p Z *-£- 


z=0 


but 


at .Li. Z 

y *.□*. = ! 

Z * X 


z=0 


therefore 



p p 

Mean square: E{x c ) - 2_ * t — 

x=0 x - 


But at x = 0 the first term ©f the series is zero, 

l ' ,r.-p. * 


therefore 


* ' ‘ »• trf 


E(X 


oo _ : ■&<: 

2) = Z 7?r4TT = 




x - l)e-V 


x=l (x " 1)1 §=1? (X “ 1)5 x=l 


— JJL X 

xe hji 
x! 


M 

7-^' r*fi‘ 


Substituting E(x) for the^iaat summation we get 

M$Wi' ...... 


-rk*. 

f-Y • "J-kf 

. -U X 


b(x*> - z ♦ St*) 

x — i mm? 1,1 


But at x * 1 the first t< 
therefore 

-u, x 
0 *11 




'i i 

JfcJ: 


the summation is zero. 


B(*2) = Z ° * n = n 2 Z n 

i=2 (x - 2)1, “ * „2 (x - 2)1 * 


x=2 

Let z = x - 2, then 




% . *3* 


‘ O* - 

E(* 2 ) = 5-St- 

z=0 *• 


+ V- 


but 


• v - 


fef# = 

z=o 


therefore 


■ • ■ • •• 
•: *»•.;•• v 


E(x 2 ) = Jl? + \L - Il({i + 1) 


2 p r ■' ' 2 

Variance: © = [\l + |i> - p, = p. 



8,5.6 Binomial Amplitude Probability Penalty Function 


Probability density function: 


PU) - inrbyT p I(1 


n 


Mean: E(x) = H x — - n - * - — P*(Lh1>): 

xl(n — x) l 

!,f 

But at x = 0 the first term of the sefieijds zero, 
therefore 

n 

z 

x=l 


f ■ 

^ -n-x 


(n - 1)1 


X^l/i : 


E(x) =np E 77’ P- ^ - P) 

c=i \ x - !)• (n - x). 


n-x 


.. «•■*! A 


Let z = x - 1 and m = n - 1, then 




E(x) = npl" £ ■■ m! p Z (l - p) m - Z 

L**o **<*- *>* & * J 


but 


^6 *'•(" - *)'• p (1 " p) 


r f 

kdS?"* 't 
#• /J8& 


therefore 

E(x) = np 
n P , 

Mean square: E(x 2 ) = 2. — p x (l - p) n " z 

x=0 x 'Cn - x)l 

But at x = 0 the first term of the serids is zero, 
therefore 7* 
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n 


SU 2 ) = Z 7 r?H 77 

x-1 (x - 1)1 (n - x)l 


p x (l - p) n ' x 


or 


_ r- (x - l)nl _x /n _,n-x 

E(x 1 - (TTijiirrioT p 11 - p) 

n . 

. v xni x,. . n-x 

+ luTT-x— p a - p) 


The second summation Is equal to E(x) = np, theref ore i \ 


n 


'|th < 

■<* 2 » - Z P x o - P)-* - ||J; 

C $ , 

" ' i 


C -1 (x - 1)1 (n - x)l 


but at x * l the f irst term of the summation Is zer 


therefore 


E(x 2 ) = n(n - l)p 2 H _\7 P X ~ 2 < 


**=2 (x - 2)1 (n - x)7 P (1 “ 


$>< '- v • *&T 


np 




Let z = x - 2 and m - n - 2, then 


•>* >• •• t'su;- 


E(x 2 ) = n(n - l)p 2 £ a .* 1 777 P*(l - P)®"* + np ¥ 

2=0 *Ma - 


but 


m 


a! 


2=0 *’•(» - *)| 


z._ .m-z , 
P (1 - P) = 1 


E(x 2 ) * n{n - l)p 2 + np 


Therefore 
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E(x 2 ) = np [(n - l )p + l] 


■ 


Variance: a 2 = np [(n - p)p + l] - (np)‘ 


a 2 = np(l - p) 



